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INTRODUCTION 

While  manual  palpation  remains  the  first  diagnostic  line  of  defense  against  breast 
cancer[l-4],  it  is  unfortunately  limited  to  relatively  large  and  superficial  lesions.  The 
central  hypothesis  of  this  work  is  that  remote  measurement  of  elasticity,  or  hardness,  of 
breast  tissues  is  possible  and  provides  unique  information  which  could  increase  detection 
and/or  characterization  of  potentially  malignant  masses  not  readily  accessible  to  manual 
exam.  Our  preliminary  studies  suggest  that  the  proposed  methods  are  capable  of 
precisely  measuring  internal  deformation  and  strain  in  three  dimensions.  These  data  are 
required  to  reconstruct  the  elasticity  distribution  within  the  object.  Consequently, 
technologies  developed  within  the  scope  of  this  project  may  have  significant  diagnostic 
value  for  the  detection  and  management  of  breast  cancer. 

Most  elastography  to  date  has  utilized  ultrasound  [5-14],  although  recently  MRI  is 
gaining  considerable  interest  [15-22].  Usually  an  external  static  or  dynamic  deformation 
is  applied  while  the  resultant  displacement  or  propagating  shear  wave  is  documented  by 
imaging  devices.  In  general,  to  reconstruct  the  tissue-specific  property  of  Young’s 
modulus  in  complex  systems  such  as  the  breast,  the  3-D  displacement  vector  must  be 
measured  over  a  3-D  volume.  In  this  project,  we  are  developing  techniques  that  measure 
the  3-D  displacement  vector  over  any  volume  in  the  object  at  high  spatial  resolution. 
These  data  are  processed  to  produce  3-D  strain  images  for  submission  to  a  3D  elasticity 
reconstruction  routine  to  map  the  relative  Young’s  modulus  within  the  volume. 

ACCOMPLISHMENTS  IN  RELATION  TO  STATEMENT  OF  WORK 

The  goal  of  this  research  program  is  to  develop  a  sensitive  diagnostic  technique 
based  on  quantitative  elasticity  imaging  for  eventual  surrogate  palpation  of  deep  lying 
breast  lesions.  To  this  end,  specific  technical  objectives/tasks  were  adopted  and 
summarized  in  the  proposal  Statement  of  Work.  Each  task  is  restated  below  along  with 
descriptions  of  relevant  Methods  and  Progress/Results. 

Technical  Objective  I:  Data  Acquisition  and  Reconstruction 

Task  1:  Month  1-3 

Construction  of  computer  controlled  hydraulic  compression  device  capable  of  producing 
an  incremental  surface  deformation  of  the  mechanical  body.  This  device  will  be 
triggered  by  NMR  imaging  system  and  have  a  simple  timing/displacement  control 
appropriate  for  phantom  studies. 

The  proposed  method  requires  that  the  imaged  object  experience  a  relatively 
minor  externally-applied  deformation  force.  Based  on  signal  optimization  simulations, 
differential  deformations  below  10%  (e.g.,  <lcm  deformation  across  a  10cm  object) 
should  be  adequate  for  generation  of  elasticity  maps.  Design  details  of  the  mechanical 
deformation  device  were  provided  in  Year  1  Progress  Report,  although  are  briefly 
summarized  here  with  the  device  schematic  is  shown  in  Figure  1.  Four  neoprene  bellows 
outside  of  the  NMR  coil  (rf-coil)  provide  both  upward  and  downward  force  to  an  acrylic 
plate  on  top  of  the  phantom.  The  top  bellows  are  on  a  common  pressure  circuit,  likewise 
the  bottom  pair  are  pressurized  together  to  yield  uniform  vertical  displacement  along  the 
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top  of  phantom.  Solid  acrylic  yolks  at  both  ends  of  the  device  insure  that  only  the  top 
plate  moves  in  response  to  pressurization.  Pneumatic  elements  are  driven  by  air-pressure 
solenoid  values  that  are,  in  turn,  controlled  by  a  transistor-transistor-logic  circuit 


Phase-encode  MRI  methods,  as  employed  here,  require  that  data  are  acquired  in 
many  segments  over  an  extended  period  of  time.  Consequently  the  mechanical 
deformation  must  be  highly  reproducible  for  each  segment.  Mechanical  stability  tests  of 
the  device  were  performed  using  a  slightly  modified  version  of  the  2D-displacement 
encoding  MRI  sequence.  The  upshot  of  these  tests  was  that  the  system  is  extremely 
stable  and  reproducible  in  application  of  a  5mm  differential  deformation  (mixing  time, 
TM=350msec).  The  observed  phase  instabilities  were  barely  above  that  allowed  by  NMR 
signal  limits.  These  results  suggest  mechanical  reproducibility  within  <20microns.  One 
significant  finding  of  these  stability  tests  was  that  the  mechanical  response  of  the 
deformation  device  was  significantly  faster  (~  150ms)  than  the  “settling  time”  of  the 
tissue-mimicking  phantoms  (~350ms).  That  is,  internal  mechanical  reflected  waves 
persist  beyond  point  the  acrylic  drive  plates  stop.  We  anticipate  that  breast  tissue  will 
have  comparably  long  settling  times,  thus  we  will  continue  to  use  long  mixing  times 
(~350ms)  in  our  experiments  despite  the  SNR  gains  of  shorter  times.  Immunity  to  long 
persistence  mechanical  waves  was  the  main  motivation  to  use  the  stimulated-echo 
approach  as  is  done  here.  In  addition,  we  propose  a  "motion-compensated"  version  of  the 
sequence  that  may  allow  use  of  shorter  mixing  times  even  during  residual  motion.  This 
concept  is  described  in  greater  detail  in  Task  11. 


Page3 


Grant  Number:  DAMD 17-97- 1-7079 


Task  2:  Month  1-6 

Development  of  various  phantoms  to  test  NMR  displacement  and  strain  data  acquisition 
and  3-D  elasticity  reconstruction  algorithms. 

Fabrication  of  phantoms  to  model  the  mechanical  and  NMR  properties  of  human 
tissue  is  an  important  step  in  development  and  verification  of  optimized  NMR 
displacement-encoded  volume  simulated-echo  pulse  sequence.  Moreover,  simple 
geometries  of  phantoms  allow  characterization  of  the  elasticity  reconstruction  algorithms. 
The  phantom  materials  should  closely  resemble  relevant  tissue  properties.  In  addition, 
these  materials  should  be  stable  (i.e.  months)  and  permit  the  non-destructive  embedding 
of  lesion-equivalent  targets  within  the  phantom. 

We  have  previously  developed  tissue-mimicking  phantoms  to  test  NMR  elasticity 
imaging  device.  These  phantoms  were  made  of  two  materials:  a)  polymer  produced  by 
M-F  Manufacturing  Co.,  Inc.  (Fort  Worth,  TX),  and  b)  gelatin  (Sigma-Aldrich  Co.,  3050 
Spruce  Street,  St.  Louis,  MO  63103  USA).  The  first  material,  plastisol,  consists  of  a 
liquid  plastic  combined  with  either  softener  (plasticizer)  or  hardener.  By  varying  the 
proportion  of  these  two  components,  it  is  possible  to  produce  composite  models  of 
desired  elasticity.  The  raw  composite  materials  were  stirred  and  heated  to  approximately 
170°C.  At  that  temperature  the  mixture  was  poured  into  molds  producing  a  tissue- 
mimicking  time-stable  phantom  of  desired  shape  and  elasticity  distribution.  However, 
two  major  deficiencies  remain.  First,  plastisol  does  not  have  desired  tissue  equivalent 
NMR  properties.  Second,  due  to  high  temperature  rise  during  phantom  preparation,  it  is 
impossible  to  produce  tissue-containing  phantoms  using  plastisol. 

Several  different  materials  were  further  considered  for  NMR  elasticity  phantoms. 
These  materials  include  Semicosil  921,  Semicosil  905  and  Silgel  612  (Wacker  Silicones 
Co.,  Adrian,  MI  49221),  and  Rhodorsil  RTV  163  (Rhone-Poulenc,  France).  From  all 
tested  materials,  the  Semocosil  921  silicone  gel  appeared  to  mimic  the  tissue  properties 
the  best.  The  elasticity  of  silicone  gel  can  be  simply  controlled  by  mixing  ratio  of  two 
components.  The  composite  material  is  then  poured  into  the  mold  and  cured  at  the  room 
temperature  enabling  fabrication  of  tissue  containing  phantoms.  These  silicone  gels  have 
shown  high  SNR  spin-echo  signal  over  the  range  of  hardness  variations.  In  addition,  the 
mechanical  properties  of  silicone  did  not  change  over  60  days,  as  was  measured  by 
Instron-type  mechanical  system.  Therefore,  it  is  possible  to  fabricate  tissue-mimicking 
phantoms  of  desired  shape  and  elasticity  distribution  using  Semicosil  921  silicone  gel. 

Based  on  favorable  NMR  and  mechanical  properties,  phantoms  were  fabricated 
using  the  Semicosil  materials.  These  phantoms  were  designed  with  simple  geometries  to 
assess  spatial  resolution  and  accuracy  of  3D  strain  imaging  3D  elasticity  reconstruction 
algorithms,  and  include  crossed  bars  of  hard  material  (6x  Young’s  modulus)  with  variable 
size  and  separation  (for  resolution)  and  simple  spherical  inclusions  for  comparison  with 
theoretical  results. 

Task  3:  Month  1-12 

Development  of  the  NMR  displacement-encoded,  volume  stimulated-echo  pulse  sequence 
optimized  for  displacement/ strain  sensitivity  and  SNR. 

Two,  3D  displacement-encode,  volume  stimulated-echo  pulse  sequences  have 
been  developed.  These  are  shown  in  Figure  2(a)  using  "classic"  3D  spatial  encoding,  and 
(b)  a  fast-spin-echo  (FSE)  extension  to  more  efficiently  encode  the  3rd  dimension.  The 
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classic  3D  sequence,  is  simplest  and  implement  but  unfortunately  it  is  clinically 
impractical  due  to  very  long  scantimes.  Incorporation  of  the  FSE  echo  train  permits  a 
scantime  reduction  by  a  factor  of  1/4,  1/8,  to  1/16  dependent  on  echo  train  length.  To 
date,  we  have  implemented  an  8-echo  FSE  sequence  where  the  FSE  phase-encoding  is 
applied  along  the  3rd  dimension  (i.e.  slice  direction).  Ideally,  all  slice  encode  steps  are 
acquired  in  one-shot  which  greatly  reduces  scan  time  (to  13min).  Unfortunately,  it 
appears  that  poor  gradient  hardware  performance  is  introducing  eddy  current  effects  that 
confound  3rd  axis  encoding  on  our  2T  MRI  system.  Phase  shifts  from  this  artifact  derail 
strain  measurement  based  on  phase.  We  are  in  the  process  of  evaluating  means  to  apply 
phase  corrections  to  reduce  imperfections  in  the  FSE  component  of  the  sequence. 
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Figure  2.  Volumetric  displacement-encoding  stimulated  echo  sequences  using  (a)  classic 
3D  encoding  and  (b)  fast-spin-echo  encoding  of  the  third  dimension. 
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Figure  3.  Volumetric  displacement-encoding  stimulated  echo  sequence  applied  to  Semicosil  rubber 
phantom  containing  crossed  ramped  bars  of  hard  material.  Shown  are  (a)  magnitude  image, 
displacement  images  in  (b)  X  right-to-left,  (c)  Y  top-to-bottom,  and  (d)  Z  through-plane.  Strain 
images  (i.e.  spatial  derivative  of  displacement)  for  X  and  Y  displacements  are  shown  in  (e)  and  (f) 
respectively. 

Images  of  tissue  mimicking  phantom  containing  ramped  bars  of  harder  material  acquired 
via  the  3D  FSE  stimulated  echo  sequence  are  shown  in  Figure  3.  Spurious  intensity 
modulations  in  the  magnitude  image  (i.e.  conventional  MR  image)  in  Figure  3(a)  indicate 
an  imperfect  reconstruction  along  the  3rd  dimension.  That  is,  there  is  an  interference  of 
signal  from  adjacent  slices.  Despite  this,  derivative  displacement  and  strain  images 
display  features  of  the  object  (Figure  3(b-f)),  however,  residual  phase  error  prevent 
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acceptable  elasticity  reconstructions.  Deformation  and  object  symmetry  meant  the 
displacement  and  strain  in  the  3rd  dimension  (ie.  perpendicular  to  the  image)  was  small. 

Computer  simulations  were  performed  to  address  the  tradeoff  between  displacement 
sensitivity  and  NMR  signal-to-noise  (SNR)  since  optimization  of  one  parameter  comes  at 
the  expense  of  another.  Key  parameters  that  affect  SNR  and  displacement  sensitivity 
include:  displacement  gradient  amplitude  and  duration;  applied  deformation  amplitude; 
image  resolution;  mixing  time,  TM;  and  tissue  water  diffusion  properties.  Results  of  the 
optimization  studies  suggest  greater  SNR  can  be  achieved  through  greater  surface  strain 
(i.e.  deformation  of  the  object  by  -10%)  with  a  corresponding  reduction  in  displacement 
encoding  (reduced  from  4  to  1  gauss/cm).  Increasing  the  applied  deformation  will  have 
the  added  benefit  to  make  the  applied  deformation  more  significant  relative  to 
background  motions  such  as  cardiac  and  respiratory.  Yet  higher  deformation,  however, 
introduce  the  need  to  consider  non-linear  elastic  effects  (see  Task  5). 


Task  4:  Month  1-12 

Expansion  of  previously  developed  linear  elasticity  reconstruction  algorithms  for 
volumetric  displacement  and  strain  NMR  measurements. 

Expansion  of  linear  elasticity  reconstruction  algorithms  were  focused  on  the 
following  two  aspects:  a)  development  of  boundary  detection  methods  to  identify  the 
regions  of  uniform  elasticity  (Young’s  modulus)  distribution,  and  b)  development  of 
volumetric  linear  elasticity  reconstruction  algorithms  for  optimized  NMR  displacement- 
encoded,  volume  simulated-echo  pulse  sequence. 

The  boundary  detection  algorithm  was  developed  for  specifically  NMR  elasticity 
imaging.  This  algorithm  is  based  on  the  stress  continuity  condition  applicable  for  soft 
tissue  deformations  used  in  these  studies.  The  relevant  boundary  definition  methods  were 
detailed  in  the  Year  1  Progress  Report. 

The  following  relates  to  extensions  to  a  3D  elasticity  reconstruction  algorithm 


(at  k<x1,x2)  (b) 


(e)  .Xj)  (d)  central  profiles 

■  1.5 

1 

0.5 
0 

-0.5 


10  20  30 
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(e)  (f)  Xj) 


<0)  Kgfx,,!^,)  (ti)  central  profiles 


10  20  30 

*  rmml 


Figure  4.  Simulated  elasticity  maps  at  planes  just  outside  of  (top  row)  and  just  within  of  (bottom 
row)  a  spherical  hard  inclusion.  The  3D  reconstruction,  k3,  (images  on  right)  was  a  more  faithful 
depiction  of  truth,  k,  (images  on  left),  relative  to  the  2D  reconstruction,  k2,  (images  in  middle). 
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wherein  complete  three-dimensional  strain  data  are  required  to  solve  for  a  general,  three- 
dimensional  object.  Previously,  only  2D  reconstruction  algorithms  have  been  applied, 
yet  these  are  anticipated  to  lead  to  errors  at  planes  where  through-plane  strain  is 
significant.  As  an  illustration.  Figure  4  shows  simulated  elasticity  distributions  at  a  slice 
just  within  (z=0.95  radius)  and  just  outside  (z=1.05  radius)  of  a  hard  spherical  inclusion. 
These  data  were  used  to  test  the  3D  elasticity  reconstruction  algorithm  and  provide 
supportive  evidence  for  the  need  of  3D  reconstruction  algorithms  in  general. 


Task  5:  Month  10-26 


Development  of  the  nonlinear  elasticity  reconstruction  model  capable  to  process  high- 
strain  NMR  images. 

Large  external  deformations  increase  the  signal  to  noise  ratio  (SNR)  of 
displacement  and  strain  images.  Unfortunately,  large  deformations  of  soft  tissue  and 
tissue-like  materials  cannot  be  described  with  a  linear  elastic  model.  A  linear  model  can 
break  down  in  two  ways.  First,  for  most  soft  tissues,  the  elastic  modulus  increases  as  a 
function  of  strain  (i.e.,  strain  hardening).  This  effect  is  often  referred  to  as  “material 
nonlinearity.”  Second,  a  more  complete  description  of  the  equilibrium  equation, 
including  non-linear  strain-displacement  relations,  must  be  used  for  large  deformations. 
This  effect  is  often  referred  to  as  “geometric  nonlinearity.”  Due  to  the  high  order 
displacement  derivatives  resulting  from  this  description,  error  propagation  must  be 
minimized  in  any  reconstruction  algorithm  using  measured  displacement  data  with  a 
finite  signal  to  noise  ratio. 

First,  we  examined  the  geometric  nonlinearities  arising  from  large  amplitude 
deformations.  This  is  demonstrated  on  the  rubber-based  phantom  with  a  single,  hard 
spherical  inclusion  (see  Figure  5).  In  general,  the  nonlinear  relation  between  the  strain 
tensor  and  the  displacement  vector,  i.e.  the  expressions  for  the  nonlinear  Lagrangian 
strain  tensor  components  are 


£ii  =  2  (M'  ’  >  +UJ  'i+ukuuk’j)>*'J  =  1-2-3 


and  do  not  depend  on  mechanical  properties  of  the  object.  In  this  equation  and  the  rest  of 
this  report,  the  lower  index  after  a  comma  means  differentiation  with  respect  to  the 
corresponding  spatial  Lagrangian  coordinate.  In  linear  case,  i.e.,  for  small  deformations, 
the  last  term  in  above  equation  is  small  and  can  be  ignored.  The  linear  and  nonlinear 
shear  components  of  the  strain  tensor  are  presented  in  Figure  5.  Note  that  these  studies 
were  conducted  on  phantom  with  almost  no  material  nonlinearities  over  the  deformation 
range  applied  (average  strain  up  to  15%).  The  specific  purpose  of  the  present  study  was 
to  explore  numerical  methods  minimizing  the  effects  of  higher  order  displacement 
derivatives  needed  to  describe  finite  amplitude  deformations  on  elasticity  reconstruction. 
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(A)  (B)  (C) 

Figure  5:  Nonlinear  shear  component  (A)  of  the  strain  tensor  is  contrasted  with  the  same 
component  computed  using  linear  model  (B),  i.e.,  under  small  deformation  assumption.  The 
difference  can  reach  almost  25%  of  the  dynamic  range  as  indicated  in  image  (C). 


Previous  algorithms  for  elasticity  reconstruction  were  formulated  using  the  set  of 
equations  describing  mechanical  equilibrium  in  a  statically  deformed,  linear  elastic 
medium.  Independent  of  the  specific  elastic  model,  however,  these  equations  can  be 
posed  in  either  differential  or  integral  form.  An  integral  representation  is  more 
appropriate  for  a  nonlinear  model  given  realistic  measurement  noise.  Numerical  methods 
have  been  developed  for  both  linear  and  nonlinear  models  exploiting  an  integral 
representation  of  the  reconstruction  equations.  Both  plane  strain  state  approach  to 
approximate  two-dimensional  displacement  and  full  three-dimensional  formulation  were 
considered. 

Second,  the  nonlinear  form  of  equilibrium  equations  using  the  nonlinear  strain- 
displacement  and  stress-strain  relationships  were  studied.  Briefly,  the  most  general 
nonlinear  mechanical  equilibrium  equations  for  a  three-dimensional  (3-D)  volume  V  of 
deformed  media  with  the  displacement  vector  U=U(x{,x2 ,x3)  =  (m1,w2,m3)  in  Cartesian 
coordinates  X=X(x,,  x2,  x})  are 

3  3 

I(Zl<V4 *=1.2,3  ■ 

j=  1  *  =  1 

Here  anj  is  a  component  of  the  2nd  ranked  stress  tensor  and  Sm  is  the  Kronecker 
delta  symbol.  This  equation  must  be  satisfied  at  every  internal  point  of  the  body.  If  the 
magnitudes  of  the  spatial  derivatives  of  all  displacement  components  are  small,  the  last 
terms  m,  ,„  can  be  omitted,  producing  the  familiar  linear  equilibrium  equations: 

3 

*'=1.2,3  . 
j- 1 

To  complete  the  system  of  equations  describing  internal  deformations,  the  relation 
between  stress  and  strain  tensors,  as  well  as  the  relation  between  the  strain  tensor  and  the 
displacement  vector,  are  needed.  Here  we  will  assume  that  the  standard  linear  relation 
between  the  stress  tensor  cr.  and  the  strain  tensor  £tj  for  incompressible  media  is  still  valid 

on  -  pSjj  +  2/ie-  , 

where  p  is  the  static,  internal  pressure  and  the  shear  elastic  modulus  fi  is  considered  a 
constant  independent  of  the  strain  magnitude.  Computing  the  spatial  distribution  of  the 
shear  elastic  modulus  is  the  goal  of  reconstruction.  Note  in  an  incompressible  material 
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such  as  soft  tissue  the  shear  and  Young’s  moduli  are  simply  proportional  (i.e.,  E  =  |i/3). 
Thus,  shear  modulus  and  Young’s  modulus  reconstructions  are  equivalent. 

For  a  general  three-dimensional  strain  state  described  by  nonlinear  equilibrium 
equation  and  stress-strain  relationship  for  incompressible  material,  the  following  system 
of  equations  can  be  obtained: 

AV(p)  =  -pB-F, 

where 


Mx,y)  = 


rl+«l,l 

M,,2 

m2’1 

1  +  u2  ,• 

VM3  »1 

«3-2 

2  ’2  M2  >3 


,  5(x,y)  =  (^)  =  V2[/, 


3  3 

F  =  (fi)  =  2{A'P  +  //[X £«U’U  +  X 0 ■ - 4 1}  > 

i'=l  i',7=l 

'P  =  (l^,)  =  (/^•).1+(Al6a),2+(^),3.  *=1.2,3; 

V  =  {d  I  dxx,d  I  dx2,d  /  (5t3),  V2  =  <72  /  <3r2  +  <?2  /  +  /  <2r2 

Note  that  the  nonlinear  form  of  the  strain  tensor  must  be  used  here. 

To  solve  this  equation  with  respect  to  unknowns  dpi  dxx,  dpi  dx2,  and  dpt  dx3 

we  first  compute  the  determinant  of  matrix  A,  which  is 


det(A)  =  1  +  DivU  +  det 


u, x  u, y 

v, x  v,y, 


+  det 


rv  v  > 

v’y  v'z 
W,y  W,zJ 


+  det 


ru,xu,z  ' 


+  det 


u, x  u,y  u,z 

v, x  v,y  v,z 

KW,X  W,yW,zJ 


where  g  =  det(g,y )  is  the  determinant  of  the  2nd  ranked  metric  tensor  gtj . 

Note  that  for  incompressible  materials  g  =  1,  anddet(A)  =  l,  which  greatly 
simplifies  the  inversion  of  (Al): 

V(p)  =  ap  +  f3, 

where  oc(x,y)  =  («;)  =  -A~lB ,  fi(x,y)  =  ($)  =  -A~lF,  /=  1,2,3- 

By  integrating  each  component  equation  of  this  equation  along  its  corresponding 
coordinate,  we  obtain: 

p(.x |  ,x2x3)  =  <Pip(Xi  ,x2x3 )  +  Fj  —  <p2p(x\ ,x2  ,  x3)  +  F2  —  (p3p(x | , x2  ,x3 )  +  F3 

where 

Xi  Xi 

<Pi  (xj  ,x2,x3)  =  exp[  J  at  (x,  ,x2,x3)dxi],  Fi(x1,x2,x3)  =  (j  Arfx,  )(pt ,  i=l  ,2,3. 

4  x?  to 

Combining  all  of  these  equations,  the  three-dimensional  case  reduces  to: 

[(^2^1,0  +0J%Lj)A>  +(<P2Fi\xo  +<P3F2\x°)  +  (F2  +  F2^  „  Pi  +2Fl  = 

+M\xo)Po  +(<PxF)\x°  +^3FiL»)  +  (F1  +F3>] 


(  <Pl  +  2F2  ~ 

mVllj  +<P2<Pl\xo)Po+(<PlF2\x°  +<P2Fl\xo)  +  (Fl+F2)l\  0<P3+2F3’ 
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where  p0  =  p(xi,x%,x%) . 

The  caveat  here  is  that  last  equation  does  not  contain  high  order  spatial 
displacement  derivatives,  compared  with  equivalent  differential  equation,  and  therefore, 
elasticity  reconstruction  should  be  more  stable.  Again,  this  equation  shows  that  spatial 
derivatives  of  all  displacement  components  are  needed  in  general  for  elasticity 
reconstruction  in  both  linear  and  nonlinear  cases,  but  any  one  of  the  three  displacement 
components  can  be  reconstructed  using  incompressibility  processing  based  on  the  relation 
(det(A)  =  1). 

Technical  Objective  II:  Phantom  Studies  on  2T  18cm  Bore  MRI  System 
Task  6:  Month  7-22 

Development  of  gel-  and  rubber-based  phantoms  with  tissue  mimicking  elastic  properties. 
Time  stable  phantoms  with  non-palpable  inclusions  of  various  shapes  and  elasticity 
contrast  positioned  at  different  locations  within  the  phantom  will  be  produced. 

During  the  first  year  of  the  project,  we  have  identified  and  tested  the  materials  for 
tissue  mimicking  phantoms.  Our  ultimate  goal  is  to  simulate  the  anatomical  and 
geometrical  features  of  the  normal  and  pathologically  transformed  breast  using  these 
materials.  The  models  of  breast  containing  single  or  multiple  inclusions  were  fabricated 
using  plastisol  material  (M-F  Manufacturing  Co.,  Fort  Worth,  Texas,  USA). 

Our  studies,  however,  suggest  that  silicone  gels  better  simulate  the  mechanical 
and  NMR  relaxation  properties  of  the  tissue.  Initially,  silicone-gel  based  homogeneous 
phantoms  were  produced  for  mechanical  and  NMR  testing  (Semicosil  921,  Wacker 
Silicones  Co.,  Adrian,  Michigan,  USA).  These  tests  showed  the  material  had  suitable 
NMR  properties  -  like  tissue.  To  control  mechanical  properties,  the  Semicosil  921 
contains  two  components,  A  and  B,  wherein  different  ratios  of  these  components  are  used 
to  vary  the  mechanical  properties  of  the  gel.  A  tissue-mimicking  phantom  was 
constructed  in  several  steps.  First,  background  material  was  prepared  by  thoroughly 
mixing  components  A  and  B  in  a  1:1  ratio,  and  then  pouring  the  mixture  into  a  154-mm 
by  80-mm  rectangular  mold.  The  mixture  was  degassed  and  cured  for  24  hours  at  room 
temperature  to  produce  a  22-mm  thick  layer.  Then  a  25-mm  diameter  hard  sphere  was 
prepared  from  a  1:2.5  mixture  of  A  and  B  and  was  placed  on  top  of  the  layer  in  the 
middle  of  the  mold.  Finally,  another  batch  of  background  material  (1:1  ratio)  was  poured 
into  the  mold  resulting  in  a  64-mm  by  80-mm  by  154-mm  phantom  with  a  single,  hard, 
spherical  inclusion  roughly  in  the  center.  At  the  same  time,  three  samples  of  each  batch 
were  taken  to  independently  assess  the  elasticity  contrast  between  the  inclusion  and 
surrounding  materials.  These  measurements  were  performed  using  the  force-deformation 
system  which  showed  that  the  inclusion  was  four  times  harder  than  the  background,  and 
that  both  background  materials  were  elastically  equivalent.  These  are  highly  stable  and 
are  still  in  use  several  months  post-fabrication. 

Task  7:  Month  12-26 

Investigation  of  the  capabilities  of  3-D  NMR  elasticity  imaging,  i.e.,  determine 
sensitivity,  accuracy  and  resolution  of  NMR  displacement  and  strain  images,  and 
reconstructed  elasticity  (Young’s  modulus)  distribution. 
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As  shown  by  simulation  presented  in  Figure  4,  elasticity  maps  reconstructed  from  3-D 
strain  data  should  be  more  accurate  than  a  2-D  elasticity  reconstructions.  This  should  be 
particularly  true  near  the  edges  of  an  object  where  the  plane-strain-state  is  violated;  that  is 
where  the  strain  along  the  “through-plane”  direction  in  non-negligible.  This  concept  was 
tested  on  silicone-gel  based  phantoms  developed  in  Task  6  using  the  volume  stimulated- 
echo  pulse  sequence  (Task  3).  A  comparison  of  the  actual  experimental  “2D” 
reconstruction  and  “3D”  reconstructed  elasticity  maps  for  2  of  32  planes  through  the 


(0  <2(x1-x2)  tel  K3^xi*x2^  (h)  Central  profiles 


10  20  30 

x,  [mm] 


Figure  6.  Reconstructed  elasticity  maps  at  planes  just  outside  of  (top  row)  and  just  within 
(bottom  row)  a  spherical  hard  inclusion  in  a  silicone-gel  phantom.  These  results  are  consistent 
with  the  simulations  (Figure  4)  that  show  the  3D  reconstruction,  k3,  (images  on  right)  more 
faithfully  depicts  the  size  of  the  inclusion. 


phantom  data  are  shown  in  Figure  6.  The  planes  selected  through  the  phantom  and 
general  format  of  Figure  6  were  chosen  for  direct  comparison  between  the  simulation 
(Figure  4)  and  experiment  (Figure  6).  Ideally,  there  would  be  no  inclusion  visible  in  the 
plane  outside  of  the  inclusion  (ie.  reconstructions  in  the  top  row  of  Figure  6  should  be 
uniform).  Violation  of  the  non-negligible  strain  remote  from  the  object,  however,  results 
in  inaccuracies  in  2D  reconstructions.  These  errors  are  clearly  reduced  in  the  3D 
reconstruction  which  confirms  simulation  results. 


Technical  Objective  III:  Translation  to  1.5T  Human  MRI  System 
Task  8:  Month  20-25 

Design  and  construction  of  the  compression  device  compatible  with  human  breast 
NMR  imaging  system  and  capable  to  produce  wide  range  of  surface  deformations. 
Development  of  sophisticated  circuits  to  control  surface  deformations  and  time 
synchronization  with  human  MRI  system. 

Translation  to  the  1.5T  unit  system  has  begun.  We  have  a  compressor  pneumatic- 
drive  system  and  constructed  the  deformation  device  suitable  for  operation  on  our 
1.5T  human  MRI  unit.  Existing  phantoms  are  used  for  testing  on  the  1.5T  unit. 
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Task  9:  Month  24-34 

Development  of  appropriate  for  clinical  studies  3-D  displacement  encoded,  volume 
stimulated-echo  pulse  sequence  on  human  MRI  system  with  data  acquisition  time 
within  or  comparable  to  regular  MRI  examination. 

Two-  and  three-dimensional  fast-spin-echo  (FSE)  are  still  being  modified  for 
displacement-sensitive,  stimulated-echo  acquisitions.  These  sequences  are  inherently 
more  complex  on  the  clinical  MRI  system,  and  as  such  are  requiring  more  time  to 
debug.  Nevertheless,  these  sequences  will  be  completed  outside  of  the  funded 
project  period. 


Task  10:  Month  22-34 

Development  of  time  efficient  elasticity  reconstruction  algorithms  more  suitable  for 
clinical  applications. 

Since  the  reconstruction  is  a  numerical  solution  to  invert  high-order  differential 
equations,  reconstruction  time  scales  with  the  acquired  quantity  of  slices  and  the  in-plane 
resolution.  Previously,  volumetric  datasets,  as  generated  in  this  project  did  not  exist,  thus 
not  all  practical  issues  related  to  3D  elasticity  reconstruction  were  known.  To  date,  most 
programming  efforts  have  been  directed  to  management  of  the  acquired  six-dimensional 
datasets  (i.e.  position  in  x,y,z,  and  displacement-sensitive  ux,uy,u2)  which  expand  to  and  a 
“12- vector”  per  pixel  (position,  and  9-elements  of  the  strain  tensor).  The  12-dimensional 
data  is  only  the  input  to  the  reconstruction. 

Both  ultrasound  and  NMR  imaging  systems  can  be  used  for  reconstructive 
elasticity  imaging.  Relative  to  ultrasound,  NMR  has  the  advantage  in  overall  resolution 
and  accuracy  for  three-dimensional  motion  estimation.  This,  in  turn,  does  not  constrain 
the  mathematical  models  used  for  elasticity  reconstruction.  Elasticity  imaging  with  NMR 
is  therefore,  volumetric  if  all  components  of  the  displacement  and  strain  are  measured. 

The  reconstruction  of  the  unknown  spatial  distribution  of  shear  elastic  modulus  in 
a  region  of  interest  (ROI)  of  three-dimensional  volume  V  in  Cartesian  coordinates  x; , 
i  =  1,2,3  under  the  assumption  of  linear  elastic  model  can  be  written  in  following  set  of 
differential  equations 

~  £ll)]>12"K/^'23)>13— ^ 

(//£j3),]|— ■(//f]3),33+[//(f33  —  )] »13-^-(^^23 )’12  ' (/^12)’23  — 

(/^23)’33"^[A(^33  —  ^22  (/^12)’13  — 

Computing  the  spatial,  3-D  distribution  of  the  shear  elastic  modulus  fi{X)  is  the  goal  of 
elasticity  reconstruction. 

Due  to  the  finite  noise  in  measured  internal  displacement  fields,  the 
computationally  more  stable  integral-based  approach  to  shear  modulus  reconstruction  is 
used,  where  corresponding  integral  equations  could  be  used  to  reconstruct  spatial 
distribution  of  unknown  shear  modulus  //(X)  =  //(xpx2,x })  instead  of  partial  differential 
equations.  For  example,  the  first  equation  can  be  written  in  integral  form  as  follows 
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S(x{yx2,x3)  —  H(Zl2  £u)  ^,l)lLo 

[M(£22  +  22  “■£ll)][co  x0  + 

J &  -  G,^.  )dx2  -  J  (G{2  -  G,'2|  )*,  =  0 , 

r°  1  r°  2 

*2  *1 

where 

<?12  =(//%), 2+(/^, 3)>3  <  G,22  =(/^12)M+(/^23),3. 

The  simplest  and  computationally  least  expensive  way  to  solve  this  set  of 
equations  is  to  reduce  any  of  those  equations  to  the  in-plane  form,  where  only  in-plane 
distribution  of  shear  or  Young’s  modulus  is  unknown  at  a  time.  It  is  possible,  for 
example,  if  for  any  region  x^<x°  the  distribution  of  shear  modulus  can  be  assumed  as 

known.  In  this  case  last  two  terms  in  the  first  equation  at  the  crossection  xi  =  y%  can  be 
replaced  by  their  finite-differences  analogs  and  final  equation  for  in-plane  shear  modulus 
distribution  /i(xu x2)  =//(*,, x2,x3)  will  be  hyperbolic.  Such  numerical  procedure  will  be 
absolutely  stable  for  any  ratio  of  grid  steps  H  t\,  ;  =  1,2  where  H  =  hi.  Therefore, 

reconstructing  volumetric  elasticity  distribution  can  be  performed  on  a  plane-by-plane 
basis.  Note  that  the  discretization  step  H  in  x3  direction  can  be  different  from  plane  to 

plane. 

A  typical  result  of  elasticity  reconstruction  is  presented  in  Figure  7.  The 
experiments  were  performed  on  the  rubber-based  block  shaped  phantom  with  a  single, 
hard  spherical  inclusion  located  approximately  in  the  middle  of  the  phantom.  In  Figure  7, 
the  elasticity  distribution  within  planes  crossing  the  center  of  the  sphere  (i.e.,  R=0)  and 
the  side  of  the  sphere  (x3  =  %  RJ  are  presented  in  Figure  8.  Also,  2-D  and  3-D  based 
reconstructions  are  contrasted.  As  expected,  the  2-D  and  3-D  reconstructions  are  very 
similar  for  the  central  plane  of  the  phantom,  where  out-of-plane  strains  are  small. 
However,  the  3-D  reconstruction  at  the  edge  of  the  inclusion  is  far  more  accurate 
compared  to  2-D  model.  The  profile  of  the  elasticity  reconstruction  throughout  the 
longitudinal  axis  is  presented  in  Figure  9,  where  again  2-D  and  3-D  results  are  contrasted. 
Clearly,  3-D  reconstruction  is  required  for  accurate  volumetric  elasticity  imaging. 
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Figure  7  (A):  2-D  elasticity  reconstruction  in  the 
central  plane  of  the  phantom  and  inclusion 


Figure  7  (B):  3-D  Elasticity  reconstruction  in  the 
central  plane  of  the  phantom  and  inclusion. 


Figure  8(A):  2-D  elasticity  reconstruction  in  the  Figure  8  (B):  3-D  elasticity  reconstruction  in  the 

plane  %  of  radius  away  from  the  inclusion  center.  plane  %  of  radius  away  from  the  inclusion  center. 
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Figure  9:  Longitudinal  elasticity  profile  of 
inclusion  using  2-D  (dotted)  and  3-D  models 
(solid).  Greater  accuracy  of  the  3-D  over  the  2-D 
reconstruction  is  apparent. 


IS  20  25  30  35 


Page  15 


Grant  Number:  DAMD 17-97-1  -7079 


Task  11:  Month  20-30 

Estimate  the  influence  of  the  cardiac  and  respiratory  motion  to  elasticity  images  and 
develop  the  approaches  to  reduce  these  artifacts. 


The  model  used  for  reconstructing  the  elastic  modulus  of  the  imaged  area  assumes  that 
the  displacement  measured  is  a  static  one.  That  is,  that  there  is  no  residual  internal 
motion  of  the  object  while  the  displacement  encoding  is  taking  place — it  is  either  in  a 
compressed  or  relaxed  state.  To  date,  all  of  our  formalism  assumed  the  object  is  static 
before  and  after  displacement  during  position  encoding  pulses.  Irreproducible 
mechanical  ring-down  of  the  deformation  system,  as  well  as  extraneous  physiologic 
motions  however,  will  cause  phase  errors.  Given  this  eventuality,  formalism  related  to 
the  displacement-encoding  was  generalized  to  include  position,  velocity,  acceleration  and 
all  other  higher-order  terms.  Mathematically,  the  encoded  phase  is  expressed  as: 

t  t  r  ~ * 

‘/>  =  rjGd(t)-r(t)  =  r]Gi(t)'  1 


r0  +  Vt  +  ■ 


Clt  + . . . 


where  phase  (</>),  the  information  accessed  with  the  imaging  sequence,  depends  on  the 
absolute  position  of  the  object  (r0),  its  velocity  (v),  its  acceleration  (a)  and  ultimately  all 
aspects  of  the  object’s  internal  motion.  So,  if  the  object  is  actually  moving  slightly 
during  the  experiment  and  we  have  data  that  contain  information  about  this  motion,  the 
following  question  arises:  if  we  use  this  data  with  a  reconstruction  that  assumes  no 
internal  motion,  what  is  the  resulting  error?  Additionally,  if  we  can  reduce  sensitivity  to 
internal  motions,  can  we  improve  the  signal-to-noise  ratio  (SNR)  of  the  map  of  the  elastic 
modulus? 


Because  the  encoding  gradient  waveform  (GJt))  is  controllable  in  the  experiments,  one  is 
able  to  reduce  or  enhance  sensitivity  to  different  aspects  of  the  imaged  object’s  motion. 
To  test  if  these  questions  would  be  of  concern  at  all,  we  developed  an  encoding  waveform 
to  provide  a  phase  image  that  is  independent  of  velocity,  but  was  otherwise  essentially 
identical  to  the  waveforms  used  in  the  past  (in  other  words,  the  phase  depended  on  all 
aspects  of  the  object’s  motion).  By  comparing  the  phase  obtained  with  the  new 
waveform  and  that  obtained  with  the  old  waveform,  one  can  see  if  the  internal  velocity  of 
the  object  contributes  to  the  encoded  phase  in  a  significant  way. 
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Velocity  compensated: 
All  motions  minus  velocity 
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Figure  10:  The  original  sequence  (on  left)  is  sensitive  to  all  motions;  the  "velocity-compensated" 
displacement  sensitive  sequence  (on  the  right)  eliminates  phase  shift  due  to  constant  velocity  motion 
during  application  of  Gd. 


Figure  1 1  shows  the  phase  images  resulting  from  both  the  old  (velocity  sensitive)  and  the 
new  (velocity  insensitive)  encoding  waveforms,  as  well  as  the  average  difference  between 
them  for  two  different  mechanical  transition  times.  The  first  transition  time  (TM  =  270ms) 
is  that  used  in  previous  experiments,  while  the  second  (TM  =  50ms)  is  a  transition  in 
which  we  know  that  the  object  was  still  moving  during  the  phase  encoding  process. 
These  results  point  to  several  things.  First,  that  the  270ms  transition  time  is  sufficiently 
long  to  allow  internal  motions  to  die  out.  Thus,  elasticity  maps  obtained  from  data  with 
this  transition  time  should  be  reasonable.  Secondly,  we  can  shorten  the  mechanical 
transition  period  at  least  until  the  two  waveforms  provide  different  results,  improving  the 
SNR  of  the  collected  data  and  the  reconstruction.  Thirdly,  it  may  be  possible  to  shorten 
this  transition  more  by  using  appropriate  encoding  waveforms,  increasing  the  SNR  and 
the  sensitivity  of  this  method  to  elasticity  variations  even  further 
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^old  ^new  A<D 


Figure  11:  Comparison  of  original  and  velocity-compensated  sequences  on  a  phantom  at  short 
mixing  time,  TM=50ms  (i.e.  where  residual  velocity  is  present),  and  long  mixing  time  TM=270ms 
(where  residual  motion  is  negligible.  The  non-zero  phase  difference  between  the  two  suggests  that 
velocity  phase  contribution  can  be  isolated,  and  thus  be  compensated. 


Task  12:  Month  30-36 

Validation  of  clinical  NMR  data  acquisition  and  elasticity  reconstruction  methods  on 
breast  tissue  mimicking  phantoms. 

This  task  is  still  incomplete.  Two-  and  three-dimensional  fast-spin-echo  (FSE)  are 
being  modified  for  displacement-sensitive,  stimulated-echo  acquisitions.  These 
sequences  are  inherently  more  complex  on  the  clinical  MRI  system,  and  as  such  are 
requiring  more  time  to  debug  than  originally  anticipated.  Nevertheless,  these 
sequences  will  be  completed  outside  of  the  funded  project  period. 


CONCLUSIONS 

Tasks  to  design,  fabricate  and  refine  the  “deformation  devices”,  essential  for 
phantom  studies,  are  complete.  The  MRI-compatible  hardware,  pneumatic  components, 
and  control  circuitry  are  now  fully  operational.  This  device  provides  excellent  control 
and  stability  in  repetitive  deformations  of  an  tissue-mimicking  objects.  Significant  effort 
was  directed  to  develop  suitable  phantom  materials.  Requirements  include:  temporally 
stable  (over  months-years);  adjustable  mechanical  properties  to  match  a  range  of  tissues; 
moldable  to  simple  geometries;  tissue-like  NMR  properties.  The  Semocosil  921  silicone 
gel  generally  meets  these  requirements  (with  the  exception  of  tissue-like  diffusion)  and 
will  be  used  in  subsequent  phantom  studies.  Three-dimensional,  stimulated-echo 
acquisition  sequences  that  have  sensitivity  to  3-dimensional  displacement  have  been 
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written  and  applied  to  gather  3D  strain  data  on  simple  objects.  These  data  have  been 
input  in  newly-designed  3D  elasticity  reconstruction  routines  to  yield,  to  the  best  of  our 
knowledge,  the  first  elasticity  reconstruction  based  on  volumetric  internal  spatial/strain 
data.  Methods  to  reduce  data  acquisition  time,  via  stimulated-echo  spatial  encoding  have 
not  yet  been  satisfactory  in  terms  of  artifact  control.  Complexity  of  fast-spin-echo 
sequences  on  the  human  MRI  unit  have  slowed  our  progress  to  translate  these  sequences 
to  the  human  system,  although  these  tasks  will  be  completed  outside  of  the  scope  of  this 
project. 
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Elasticity  Reconstructive  Imaging  by  Means  of 
Stimulated  Echo  MRI 

Thomas  L.  Chenevert,  Andrei  R.  Skovoroda,  Matthew  O’Donnell, 
Stanislav  Y.  Emelianov 


A  method  is  introduced  to  measure  internal  mechanical  dis¬ 
placement  and  strain  by  means  of  MRI.  Such  measurements 
are  needed  to  reconstruct  an  image  of  the  elastic  Young’s 
modulus.  A  stimulated  echo  acquisition  sequence  with  addi¬ 
tional  gradient  pulses  encodes  internal  displacements  in  re¬ 
sponse  to  an  externally  applied  differential  deformation.  The 
sequence  provides  an  accurate  measure  of  static  displace¬ 
ment  by  limiting  the  mechanical  transitions  to  the  mixing 
period  of  the  simulated  echo.  Elasticity  reconstruction  in¬ 
volves  definition  of  a  region  of  interest  having  uniform 
Young’s  modulus  along  its  boundary  and  subsequent  solution 
of  the  discretized  elasticity  equilibrium  equations.  Data  acqui¬ 
sition  and  reconstruction  were  performed  on  a  urethane  rub¬ 
ber  phantom  of  known  elastic  properties  and  an  ex  vivo  ca¬ 
nine  kidney  phantom  using  <2%  differential  deformation. 
Regional  elastic  properties  are  well  represented  on  Young’s 
modulus  images.  The  long-term  objective  of  this  work  is  to 
provide  a  means  for  remote  palpation  and  elasticity  quantita¬ 
tion  in  deep  tissues  otherwise  inaccessible  to  manual  palpa¬ 
tion. 

Key  words:  elastic  Young’s  modulus;  magnetic  resonance 
imaging;  elastography;  strain  imaging. 

INTRODUCTION 

It  is  well  known  that  tissue  elastic  properties  may  be 
altered  by  tumors.  Young’s  elastic  moduli  may  differ  by 
orders  of  magnitude  in  soft  tissues  in  various  physiologic 
states  (1,  2).  This  finding  is  the  physical  basis  behind 
manual  palpation  used  to  detect  “hard”  masses  (3,  4). 
Indeed,  physical  examination  is  the  first  diagnostic  line 
of  defense  against  breast  cancer,  because  nodule  hard¬ 
ness  raises  suspicion  of  malignancy.  Detection  of  a  new 
breast  mass  by  physical  examination  is  often  sufficient 
for  surgical  excisional  biopsy,  even  when  not  corrobo¬ 
rated  by  other  diagnostic  tests.  Manual  palpation  of  the 
prostate,  superficial  lymph  nodes,  and  abdominal  organs 
are  also  commonly  performed.  Unfortunately,  sensitivity 
of  palpation  is  relatively  poor  within  deep,  dense,  or 
heterogeneous  tissues.  Although  the  touch  of  a  skilled 
interpreter  is  considered  a  powerful  diagnostic  instru¬ 
ment,  most  lesions  detected  by  palpation  tend  to  be  rel¬ 
atively  large  and  superficial. 
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Scientists  are  attempting  to  electronically  extend  the 
touch  of  the  physical  examiner  using  a  variety  of  image- 
based  techniques  that  infer  tissue  elasticity.  The  essential 
element  is  measurement  of  internal  motion  and  strain  in 
tissue  structures  experiencing  mechanical  stress.  To 
date,  most  “elastography”  has  used  ultrasound  to  track 
the  relative  motion  of  targets  by  specular  reflection  (5-7), 
by  Doppler  techniques  (8-10),  by  cross-correlation  of 
raw  or  processed  acoustic  echoes  (11,  12),  or  by  tracking 
speckle  patterns  (13-15).  Usually  an  external  static  or 
dynamic  deformation  is  applied  while  internal  displace¬ 
ments  or  propagating  shear  waves  are  documented  by 
imaging. 

MRI  has  also  been  used  to  measure  internal  displace¬ 
ment  and  strain  components  of  the  heart  using  spatial 
magnetization  tagging  (16,  17)  and  phase-based  velocity 
encoding  (18).  Elasticity  reconstruction  of  an  externally 
deformed  phantom  was  demonstrated  using  magnetiza¬ 
tion  tagging,  but  this  method  has  spatial  resolution  lim¬ 
ited  by  the  tagged  grid  size  and  only  measures  2D  motion 
(19).  More  recently,  motion  phase  encoding  by  means  of 
bipolar  gradients  was  used  to  produce  two-dimensional 
(2D)  displacement  and  strain  maps  in  media  mechani¬ 
cally  driven  by  external  forces  (20-23).  Strain  and  dis¬ 
placement  maps  infer  internal  elasticity  but  are  also 
strongly  affected  by  the  applied  deformational  geometry. 
Consequently,  these  maps  do  not  uniquely  reflect  inter¬ 
nal  tissue  properties  (i.e.,  elastic  Young’s  modulus). 
Maps  of  dynamic  strain-wave  propagation,  however,  do 
allow  measurement  of  local  strain  wavelength  or  velocity 
from  which  the  local  elastic  modulus  can  be  derived  (20). 
Shear-wave  attenuation,  interference  horn  standing 
waves  off  multiple  reflectors,  and  limited  resolvable 
points  over  the  shear  wavelength  are  potential  drawbacks 
of  this  approach. 

Relative  to  ultrasound,  MRI  has  an  advantage  in  overall 
resolution  and  accuracy  for  multidimensional  displace¬ 
ment  and  strain  measurement  needed  for  elasticity  re¬ 
construction.  Ultrasound  can  accurately  measure  axial 
(i.e.,  along  the  beam  axis)  motion  at  high  spatial  resolu¬ 
tion  (< millimeter),  but  lateral  displacement  is  measured 
at  much  lower  spatial  resolution  defined  by  the  depth- 
dependent  acoustic  beam  width.  The  third  dimension  is 
generally  not  even  considered  given  the  limitations  of 
ultrasound.  Consequently,  reduced  motion  dimensional¬ 
ity  and  overall  low  motion  resolution  of  the  imaging 
system  compromise  the  “elastogram.”  This  shortcoming, 
in  turn,  constrains  the  mechanical  model  used  in  elastic¬ 
ity  reconstruction.  To  date,  only  ID  motion  models  have 
been  applied  to  ultrasound-derived  elasticity  images  of 
tissues  in  vivo  (12,  24).  More  accurate  elasticity  images 
are  achieved  by  properly  controlling  external  deforma- 
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tions,  leading  to  2D  elasticity  reconstruction  within  the 
imaging  plane  (25,  26). 

In  this  work  we  present  a  method  to  spatially  encode 
internal  displacement  of  an  object  that  has  undergone  an 
externally  applied  “static”  deformation  with  subsequent 
reconstruction  into  elasticity  maps.  Unlike  dynamic 
techniques  designed  to  estimate  elasticity  from  observa¬ 
tions  related  to  strain-wave  propagation,  static  elasticity 
reconstruction  involves  estimation  of  local  strain  from 
displacement  and  numerical  solution  of  differential  elas¬ 
ticity  equilibrium  equations. 


THEORY  OF  RECONSTRUCTIVE  ELASTICITY 
IMAGING 


The  goal  of  elasticity  imaging  is  to  reconstruct  the  elastic 
modulus  of  a  desired  tissue  region  using  available  mea¬ 
surements  of  displacement  and  strain  components.  In¬ 
deed,  the  mechanical  properties  of  tissue  are  ultimately 
linked  to  the  patterns  of  internal  deformations,  but  the 
deformational  geometry  can  greatly  affect  these  patterns 
as  well.  To  uniquely  image  tissue  elasticity,  the  Young’s 
modulus  must  be  reconstructed  from  estimates  of  inter¬ 
nal  displacement  and  strain. 

In  this  paper,  the  general  approach  to  elasticity  recon¬ 
struction  was  based  on  a  model  of  linear,  elastic,  isotro¬ 
pic,  incompressible  media  (26,  27).  The  key  equations 
and  considerations  are  briefly  presented  here.  A  more 
detailed  description  of  elasticity  reconstruction  is  given 
in  an  earlier  publication  (26). 

In  linear  elasticity,  the  components  of  the  strain  (e^) 
and  stress  (o^)  tensors  under  static  deformation  are: 


1  (  dll;  dizA 

2  \dXj  dXj) 


1 

=  -(uifJ  +  u}ii) 


[1] 
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<Tj}  =  p8jj  H-  -  Eejj  [2] 

where  u{  is  a  component  of  the  displacement  vector  U  = 
(ua,  u2,  u3)  in  Cartesian  coordinates  X  =  [xlf  x2,  x3),  p  is 
the  static  internal  pressure,  8^  is  the  Kronecker  delta 
symbol,  and  E  =  E[xlf  x2 ,  x3 )  is  the  Young’s  elastic 
modulus.  Note  in  Eq.  [1],  and  the  entire  paper,  the  lower 
index  after  a  comma  means  differentiation  with  respect 
to  the  corresponding  spatial  coordinate. 

The  static  deformation  of  the  medium  can  be  described 
by  the  equilibrium  condition: 


3 

2  <r,N  +  //  =  0,  2  =  1,2,3  [3] 

where  ft  is  the  body  force  per  unit  volume  acting  in  the  x{ 
direction.  In  addition,  volume  conservation  for  an  in¬ 
compressible  medium  leads  to  the  following  relationship 
between  displacement  and  strain  components 

V  •  U  =  Sn  +  £22  “h  £33  =  LZlJl  ^2>2  "h  ^3>3  =  0  [4] 

Using  Eqs.  [l]  and  [2]  for  stress  and  strain  components, 
and  the  incompressibility  Eq.  [4],  the  equilibrium  condi¬ 


tion  with  eliminated  internal  pressure  p  can  be  rewritten 
in  the  following  form: 

2s12(E,11  —  E, 22)  +  2(u2>2  “  ^i>i) E,i2  “h  2£23E,13 
—  2e13E,23  +  (V2U2  +  <*>12>l)  fi>l  “  (V2Ua  —  0)12,2)E,2 
+  0)12,3E,3  4-  V2G)12E  +  3 (/25l  “  ft, 2)  =  0 


2eia(*,u  E,33)  +  2 £23-E ,12  4-  2(113,3  ^i»i)^»i3 

—  2s12E,23  +  (V2u3  4-  co13,1)£’,1  +  a)13,2E,2  [5] 

“  (V2Ui  —  co13,3)E,3  4-  V2co13E  4-  3 (/3)1  _/i ,3)  =  0 

2623(^22  —  £>33)  +  2s13E,12  —  2 e12E,13  4-  2 ( lz3 , 3 

U2» 2)^523  "h  (W23j1-^'>1  “f  (V2U3  +  (023,2)E,2  ~  (V2U2 

~~  (*)23>3)E>3  ’h  ^2<y23-^'  "f  3(/3,2  “  /As)  “  0 

m12  =  ^2»1  “  ^1>2»  w13  =  U3>1  “  ul>3>  ^23  =  U3>2  “  U2>3 

Clearly,  the  elasticity  reconstruction  process  based  on 
Eq.  [5]  requires  accurate  measurements  of  the  displace¬ 
ment  vector,  or,  to  be  more  precise,  requires  accurate 
estimation  of  up  to  third-order  spatial  derivatives  of  the 
displacement.  Equation  [5]  can  also  be  written  in  terms  of 
spatial  derivatives  of  strain  tensor  components.  The 
unique  solution  of  the  system  of  coupled  partial  differ¬ 
ential  equations  [5]  is  determined  by  the  boundary  con¬ 
ditions,  i.e.,  the  elastic  modulus  E[X)  must  be  given  along 
the  boundary  of  the  reconstruction  region  of  interest 
(ROI).  It  should  also  be  noted  here  that  the  analytical 
solution  of  Eq.  [5]  is  not  generally  possible,  and  numer¬ 
ical  methods  must  be  developed  to  solve  this  system  of 
partial  differential  equations. 

Based  on  the  particular  geometry  of  the  phantoms  and 
deformation  system  used  in  these  experiments,  Eq.  [5] 
was  simplified.  A  2D  approximation  of  Eq.  [5]  was  used 
because  the  imaged  plane  was  near  the  center  of  the 


FIG.  1.  Stimulated-echo  data  acquisition  and  object  deformation 
sequence.  Mechanical  transitions  occur  during  the  long  mixing 
time  (TM)  such  that  a  static  displacement  equilibrium  is  achieved. 
Local  displacement  between  deformation  states  “A”  and  “B”  are 
encoded  by  phase  shift  proportional  to  “G-displace”  amplitude 
Gd,  and  duration  r. 
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FIG.  2.  Urethane  rubber  phantom  containing 
two  hard  inclusions  and  held  between  defor¬ 
mation  plates  on  top  and  bottom  of  the  phan¬ 
tom.  Images  of  stimulated  echo  (a)  magnitude, 
and  phase  shift  proportional  to  (b)  vertical  and 
(c)  horizontal  differential  displacement.  Dis¬ 
placement  sensitivity  of  15.33  7r/mm  and  dif¬ 
ferential  displacement  of  «1.4  mm  between 
deformation  plates  was  used. 
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phantom.  For  such  a  plane  (arbitrarily  denoted  as  x3  =  0), 
the  displacement  vector  components  do  not  vary  signif¬ 
icantly  as  a  function  of  the  out-of-plane  x3  coordinate, 
and  therefore,  the  “plane  strain  state”  condition  is  appli¬ 
cable.  With  this  condition,  Eq.  [5]  reduces  to  a  single 
nontrivial  equation: 


founded  by  interference  of  the 
primary  shear  wave  with  re¬ 
flected  or  standing  shear 
waves.  To  avoid  this  condi¬ 
tion,  a  “static”  displacement 
encoding  approach  was 
adopted.  It  requires  measure¬ 
ment  of  internal  displacement 
between  two  or  multiple  de¬ 
formations  while  the  object  is 
in  mechanical  equilibrium  for 
each  measurement.  A  stimu¬ 
lated  echo  sequence  with  dis¬ 
placement-encoding  gradient 
pulses  is  used  to  achieve  this, 
as  shown  in  Fig.  1.  Mechanical 
transition  from  state  “A”  to 
state  “B”  occurs  during  the 
stimulated  echo  mixing  time, 
TM.  A  relatively  long  mixing 
time  allows  long-lived  elastic 
vibrations  to  dampen  before 
spatial  encoding.  Because  the 
relevant  magnetization  during 
TM  is  longitudinal,  it  is  unaf¬ 
fected  by  potentially  ill-de¬ 
fined  motions  during  the  me¬ 
chanical  transition  period.  As 
a  result,  a  more  accurate  static 
deformation  measurement  is 
achieved.  Also  note  that  pre¬ 
cise  synchronization  of  me¬ 
chanical  and  pulsed  gradient 
events  is  not  critical  as  long  as 
the  mechanical  transition  be¬ 
gins  after  the  second  RF  pulse 
and  is  complete  before  the  third  RF  pulse.  Similarly,  a 
long  delay  in  TE  could  be  used,  but  this  is  done  at  the 
expense  of  signal  lost  to  T2  decay. 

Local  displacement  is  encoded  by  means  of  phase  shift 
governed  by  pulsed-field  gradient  factors, 

=  y  Gd  T  [7] 


(£>n  -  £>22)(Ul>2  +  U2,l)  +  2£>12(u2,2  -  Ulfl) 

+  2 Efl(u2,n  4-  u 2 , 2 2 )  -  2 E,2(uull  +  Uu 22)  ^ 

+  E(u2, 111  +  U2, 221  —  ^1  >112  ~  Ul>222)  d"  3(/2>1  ~/i,2)  =  0 

Under  conditions  in  which  these  assumptions  are  valid, 
only  in-plane  displacement  ua  and  u2  or  their  derivatives 
are  needed  to  reconstruct  the  modulus  in  the  plane  x3  =  0. 

Static  Displacement  Measurement  by  Means  of 
Stimulated  Echo  MRI 

Shear-wave  propagation  speed  in  soft  tissue  is  1-20  m/s. 
Consequently,  a  shear  wave  imparted  by  a  single-stroke 
or  an  oscillating  deformation  force  may  require  tens  of 
milliseconds  to  traverse  an  object  ^100  mm  in  size.  The 
time  for  reflected  waves  to  dampen  may  be  much  longer. 
“Dynamic”  measurements,  which  encode  displacement 
during  shear-wave  propagation,  are  potentially  con¬ 


where  displacement  sensitivity,  <bd,  is  in  units  of  (radi¬ 
ans/distance).  A  phase  reference  acquisition  is  required 
for  each  displacement  encode  condition  to  remove  pre¬ 
existing  phase  shifts  that  are  unrelated  to  displacement. 
Reference  data  are  acquired  using  the  same  pulse  se¬ 
quence,  including  displacement  encode  gradient  pulses, 
but  with  the  object  maintained  in  state  B.  Note  that  all 
spatial  encoding  occurs  from  the  third  RF  pulse  and 
beyond.  That  is,  the  object  is  in  one  deformation  state 
(state  B)  for  both  displacement  and  phase  reference  ac¬ 
quisitions  during  all  spatial  encoding  segments  of  the 
sequence.  Consequently,  image  registration  or  feature 
tracking  algorithms  (14)  are  not  required  to  estimate  dis¬ 
placement.  Instead,  local  displacement  is  encoded  di¬ 
rectly  by  local  phase  of  a  corrected  dataset,  Scor,  given  by, 


SCor(r) 
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FIG.  3.  Strain  images  calculated  from  first-order  derivatives  of  the  displacement  images  repre¬ 
sented  in  Fig.  2.  Normal  strains  (a)  and  (b)  e22,  and  (c)  shear  strain  e12  =  e21  reflect  internal 
elastic  properties  and  the  externally  applied  deformation  field.  The  plane  strain  state  assumption 
and  phantom  incompressibility  suggest  e-,-,  ~  ~e22  as  is  supported  by  the  relatively  featureless 
map  of  su  +  e22  in  (d). 


where  SA  and  SB  are  the  data  acquired  with  the  object 
initially  in  state  A  and  state  B,  respectively.  The  un¬ 
wrapped  phase  of  the  corrected  dataset  and  Eq.  [7]  pro¬ 
vides  a  local  measure  of  displacement,  U,  by  means  of; 

<p(f)  =  <bd  •  [rA  -  rB]  =  $>d  *  U(r)  [9] 

Most  sources  of  phase  error,  such  as  static  field  inho¬ 
mogeneity,  tend  to  be  slowly  varying  functions  of  posi¬ 
tion.  Therefore,  the  phase  reference  datasets  may  be  ac¬ 
quired  at  relatively  low  spatial  resolution  to  reduce  scan 
time. 

METHODS 

Data  Acquisition 

Elasticity  imaging  was  performed  on  two  phantoms.  One 
was  an  85-mm  diameter  cylindrical  urethane  rubber 
phantom  containing  two  8-mm  cylinders  of  hard  mate¬ 
rial.  Previously,  the  ratio  of  Young’s  modulus  between 
the  inclusion  and  background  material  was  measured 


10.5  ±  1.5  for  30%  surface  de¬ 
formation  (19).  In  the  present 
study,  smaller  surface  defor¬ 
mation  was  applied,  and  con¬ 
sequently,  the  elasticity  con¬ 
trast  within  this  phantom 
should  be  less  than  the  previ¬ 
ously  measured  ratio.  Whereas 
the  mechanical  properties  of 
the  rubber  phantom  mimic 
soft  tissue,  the  phantom  had 
inherently  poor  NMR  signal.  A 
more  tissue-equivalent  phan¬ 
tom  in  terms  of  NMR  and  me¬ 
chanical  properties  was 
achieved  by  embedding  a  fresh 
canine  dog  kidney  (<24  h  ex 
vivo)  into  a  130  mm  X  105 
mm  X  75  mm  block  of  5%  gel¬ 
atin.  One  hour  before  MRI,  10 
ml  of  5%  glutaraldehyde  solu¬ 
tion  was  injected  into  kidney 
parenchyma  to  create  a  hard 
lesion(s). 

Phantoms  were  held  se¬ 
curely  in  place  under  moder¬ 
ate  preload  pressure  by  two 
parallel  acrylic  plates.  When 
pressure  to  the  top  plate  was 
released,  the  phantom  recoiled 
vertically  an  amount  con¬ 
strained  by  physical  stops. 
Maximum  vertical  displace¬ 
ment  was  <1.5  mm,  which 
represented  <2%  differential 
between  state  “A”  (greater  de¬ 
formation)  and  state  “B”  (less 
deformation) .  D  eformation 

was  actuated  pneumatically 
by  an  air-filled  bladder  on  top 
of  the  phantom  holder.  Pneu¬ 
matic  pressure  was  stepped  by  a  remote  solenoid  valve 
with  timing  controlled  by  an  external  transistor  transistor 
logic  (TTL)  gate  circuit  triggered  by  the  pulse  sequence. 

Displacement  encoding  gradient  pulse  duration,  t  = 
4.5  msec,  and  amplitude,  Gd  =  40  mT/m,  provided  a 
displacement  sensitivity  of  <hd  -  15.33  7r/mm  by  means 
of  Eq.  [7].  The  displacement  encoding  direction  was  al¬ 
ternated  each  pulse  repetition  between  vertical  and  hor¬ 
izontal.  For  the  urethane  rubber  phantom,  acquisition 
parameters  were  TR  -  1.3  s,  TM  =  350  msec,  TE  =  50 
msec,  128  X  128  matrix,  four  signal  averages,  100  mm 
field  of  view,  and  10-mm  section  thickness.  An  addi¬ 
tional  128  X  32  matrix  acquisition  was  collected  for 
phase-reference  correction  of  vertical  and  horizontal  en¬ 
coded  data.  The  kidney  phantom  acquisition  parameters 
were  TR  =  Is,  TM  =  200  msec,  TE  —  76  msec,  256  X  256 
matrix,  two  signal  averages,  150  mm  field  of  view,  and 
5-mm  section  thickness;  with  a  256  X  32  dataset  acquired 
for  phase  correction.  All  experiments  were  performed  on 
a  2  T  18 -cm  bore  MRI  system  (Bruker,  formerly  GE  NMR 


*486 


Chenevert  et  al. 


ary  conditions.  Therefore,  a  rectangular  ROI  was  identi¬ 
fied  within  the  imaging  planes  for  both  phantoms.  For 
the  phantom  with  two  hard  inclusions,  the  ROI  was  a 
region  of  63  X  50  mm  positioned  approximately  in  the 
center  of  the  phantom  and  included  both  inclusions.  For 
the  canine  kidney  phantom,  the  rectangular  94  X  51  mm 
ROI  included  the  whole  kidney  cross-section.  In  both 
cases,  the  Young’s  modulus  value  along  the  ROI  boundary 
was  set  to  “one”  resulting  in  reconstruction  of  relative 
Young’s  modulus.  More  detailed  analysis  and  discussion  of 
defining  the  ROI  was  considered  previously  (26). 

Elasticity  reconstruction  Eqs.  [5]  and  [6]  assume  that 
spatial  derivatives  of  the  Young’s  modulus  are  continu¬ 
ous  functions.  To  ensure  continuous  elasticity  distribu¬ 
tion,  the  spatial  derivatives  of  the  displacement  were 
low-pass  filtered  before  elasticity  reconstruction,  result¬ 
ing  in  mild  spatial  resolution  reduction. 

After  defining  the  boundary  conditions,  Eq.  [6]  was 
discretized  over  the  ROI  with  the  same  grid  spacing  as 
the  MR  images,  where  all  spatial  derivatives  of  the  dis- 
placement/strain  (i.e.,  coefficients  in  Eq.  [6]  for  unknown 
Young’s  modulus  distribution)  were  approximated  „ by 
finite  differences.  The  linear  set  of  equations  resulting 
from  discretization  of  Eq.  [6]  was  solved  iteratively, 
where  the  error  in  each  step  was  estimated  by  averaging 
the  left-hand  side  of  Eq.  [6]  over  the  ROI  using  the  current 
estimate  of  the  elasticity  distribution.  From  step  to  step, 
the  Young’s  modulus  distribution  was  updated  based  on 
the  changes  in  the  average  error. 

RESULTS 


FIG.  4.  (a)  Reconstructed  map  of  Young’s  modulus  for  the  urethane 
rubber  phantom  within  a  63  x  50  mm  region.  The  boundary  of  the 
region  is  defined  to  have  a  Young’s  modulus  =  1 .  Relative  Young’s 
moduli  for  inclusions  and  background  material  are  shown  in  (b). 

Instruments),  using  a  150-mm  transmit/receive  birdcage 
coil. 

Data  Processing 

Time-domain  data  were  transferred  for  off-line  process¬ 
ing  as  follows.  Phase  reference  datasets  were  zero-filled 
and  2D  Fourier  transformed  to  a  128  X  128  or  256  X  256 
matrix  for  phase  correction  by  means  of  Eq.  [7].  The 
resulting  phase  maps  were  used  to  estimate  the  spatial 
derivatives  of  the  in-plane  displacements  necessary  for 
elasticity  reconstruction  (6).  Note  phase  unwrapping  is 
not  strictly  required  because  only  phase  derivatives  are 
used.  Assuming  the  displacement  fields  are  continuous, 
resulting  in  small  differential  displacement  at  any  pixel 
compared  with  the  total  displacement,  the  differential 
displacement  between  two  neighboring  pixels  was  di¬ 
rectly  computed  from  the  angle  of  the  complex  multipli¬ 
cation  of  each  pixel  with  the  conjugate  of  the  neighboring 
pixel,  then  scaled  by  lM>d. 

Solving  Eq.  [6]  for  unknown  E  (xl,  x2)  performed  the 
elasticity  reconstruction,  i.e.,  reconstruction  of  the  spa¬ 
tial  distribution  of  elastic  Young’s  modulus.  As  was 
noted  previously,  the  unique  solution  of  a  boundary 
value  problem  (Eq.  [5]  or  [6])  is  determined  by  the  bound¬ 


Magnitude  and  corrected  phase  images  of  the  urethane 
rubber  phantom  are  shown  in  Fig.  2.  Given  that  <f>d  = 
15.33  7r/mm,  the  number  of  2 tt  phase  bands  in  Fig.  2b 
indicates  that  the  vertical  excursion  of  the  phantom  top 
relative  to  the  bottom  was  ^1.4  mm;  similarly,  the  rela¬ 
tive  lateral  displacement  between  left  and  right  edges  of 
the  phantom  was  ~0.8  mm  (Fig.  2c).  Reduced  phase 
slopes  in  the  regions  of  the  hard  inclusions  are  clearly 
visible  on  the  phase  images.  Normal  strain,  e11  and  e22, 
and  shear  strain,  e12  =  e21,  maps  are  illustrated  in  Figs. 
3a— 3c,  respectively.  The  observed  contrast  reversal  be¬ 
tween  £ia  and  e22  is  a  result  of  phantom  incompressibil¬ 
ity  (like  soft  tissue),  which  yields  =  -e22  assuming 
negligible  out-of-plane  strain.  Consequently,  the  sum  (elt 
+  e22)  is  relatively  “flat”  as  shown  at  equivalent  gray¬ 
scale  settings  in  Fig.  3d.  Also  note,  whereas  the  strain 
maps  clearly  exhibit  object-specific  detail  (i.e.,  inclu¬ 
sions),  features  related  to  the  applied  external  deforma¬ 
tion  are  quite  conspicuous.  This  fact  demonstrates  why 
an  elasticity  reconstruction  is  needed.  The  Young’s  mod¬ 
ulus  was  reconstructed  for  a  63  X  50  mm  region  as 
presented  in  Fig.  4a.  The  boundary  of  the  elasticity  re¬ 
construction  area  was  defined  to  have  value  “1.”  The 
relative  elastic  moduli  for  select  regions  are  indicated  in 
Fig.  4b  and  are  consistent  with  the  known  elasticity  of 
these  materials  (19). 

The  canine  kidney  phantom  is  shown  in  Fig.  5  as 
magnitude  (Fig.  5a),  vertical  (Fig.  5b),  and  horizontal 
(Fig.  5c)  phase  shift  images.  It  is  apparent  from  the  ver¬ 
tical  phase  image  (Fig.  5b)  that  although  ^1.24-mm  rel- 
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FIG.  5.  Stimulated  echo  images  of  canine  kid¬ 
ney  imbedded  in  gelatin,  (a)  Magnitude  image, 
and  (b)  vertical  and  (c)  horizontal  displacement 
phase  images  in  which  <1.3  mm  of  vertical 
differential  deformation  and  15.33  7r/mm  dis¬ 
placement  sensitivity  were  applied. 
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ative  displacement  spans  the  full  phantom,  only  ~0.26 
mm  of  it  is  within  the  kidney.  Consequently,  there  is  a 
high  concentration  of  strain  near  the  gel-kidney  interface, 
as  is  clear  on  the  strain  maps  (Fig.  6).  As  noted  before, 
strain  images  exhibit  contrast  related  to  a  combination  of 
internal  structure  and  the  externally  applied  deforma¬ 
tion.  Elasticity  reconstruction,  however,  reduces  the  am¬ 
biguity  and  exhibits  contrast  dominated  by  internal  elas¬ 
tic  properties  (Fig.  7a).  It  is  also  encouraging  to  note  that 
whereas  there  was  only  moderate  strain  contrast  within 
the  kidney,  elasticity  contrast  within  renal  parenchyma 
and  central  sinus  are  well  distinguished  on  the  Young’s 
modulus  image.  Moreover,  the  site  of  glutaraldehyde  in¬ 
jection  (top-right  quadrant  of  images)  exhibited  the  high¬ 
est  relative  Young’s  modulus.  Approximately  20  h  after 
MRI,  the  kidney  phantom  was  sliced  at  a  plane  corre¬ 
sponding  to  that  studied  by  MRI.  An  optical  image  of  this 
slice  is  shown  in  Fig.  7c.  The  freshly  cut  surface  was  pal¬ 
pated  such  that  areas  of  relatively  “hard”  parenchyma 
could  be  noted.  Arrows  in  Fig.  7c  mark  the  most  conspic¬ 
uous  areas  of  hardness;  the  largest  area  corresponds  to  the 
high  Young’s  modulus  region  in  the  upper-right  quadrant 
of  the  kidney. 

DISCUSSION 

In  this  work  we  introduce  a  method  to  image  and  quan¬ 
tify  internal  elastic  properties  of  an  object  by  means  of 
displacement-sensitive  MRI  with  associated  elasticity  re¬ 
construction.  The  data  acquisition  segment  employs  gra¬ 
dient  pulses  to  encode  internal  displacement  by  means  of 


phase  using  a  stimulated  echo 
sequence.  Internal  displace¬ 
ments  occur  in  response  to  an 
external  deformation  force 
synchronized  to  the  acquisi¬ 
tion  sequence.  By  timing  con¬ 
trol,  mechanical  motion  oc¬ 
curs  while  the  relevant 
magnetization  is  longitudinal. 
The  stimulated  echo  allows 
extension  of  the  mechanical 
transition  period  to  avoid  po¬ 
tentially  long-lived  or  ill-de¬ 
fined  oscillations  within  the 
object  such  that  an  estimate  of 
“static”  displacement  is 
achieved.  Image  registration  or 
feature  tracking  (14)  is  not  re¬ 
quired  because  the  object  is  in 
one  deformation  state  for  all 
spatial  encoding.  Tx  relaxation 
and  diffusion,  which  erode 
signal,  ultimately  set  practical 
limits  on  this  period.  In  these 
experiments,  200  to  350  msec 
was  sufficient  to  allow  static 
displacement  measurement  of 
rubber  and  gelatin/tissue 
phantoms  using  a  simple  air- 
bladder  pneumatic  system. 
This  arrangement  required  the 
object  to  “passively”  recoil  during  the  TM  period. 
Clearly,  a  faster  deformation  system  can  be  built  using 
external  forces  to  “actively”  deform  the  object  during  the 
stimulated  echo  mixing  period.  A  shorter  mixing  time 
would  then  be  used  yielding  a  larger  signal. 

To  date,  two  approaches  are  present  in  elasticity  im¬ 
aging:  static  reconstructive  elasticity  imaging  (25,  26)  and 
dynamic  shear-wave  elasticity  imaging  (20,  21,  23).  In 
both,  an  external  static  or  dynamic  deformation  is  ap¬ 
plied  while  the  resulting  displacement/strain  or  propa¬ 
gating  shear  wave  is  detected  using  an  imaging  modality. 
In  reconstructive  elasticity  imaging,  the  elasticity  distri¬ 
bution  must  be  reconstructed  from  static  displacement 
and  strain  images.  The  ability  to  control  the  internal 
deformation  pattern  by  varying  the  externally  applied 
load,  and  high  SNR  displacement  and  strain  estimates 
are  the  benefits  of  this  method,  although  numerical  re¬ 
construction  algorithms  are  required.  Wherein  the  mod¬ 
eled  assumptions  are  valid,  these  algorithms  exist.  For 
more  general  applications,  they  must  be  refined  further. 
In  shear-wave  elasticity  imaging,  local  shear  wavelength 
measurements  allow  direct  and  simple  calculation  of  the 
shear  elastic  modulus.  However,  the  interference  of  shear 
waves  reflected  from  any  elasticity  inhomogeneities 
within  the  tissue,  along  with  attenuation  of  shear  waves, 
and  conversion  between  shear  and  bulk  waves  are  chal¬ 
lenges  of  this  method. 

In  these  experiments  a  displacement  sensitivity  of  <£>d 
=  15.33  7r/mm  was  achieved  using  moderate  gradient 
factors.  The  ultimate  quality  of  Young’s  modulus  recon¬ 
struction  depends  on  the  induced  phase  shift,  equal  to 
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FIG.  6.  Strain  maps  of  the  kidney  phantom. 
Normal  strains  (a)  eu  and  (b)  e22,  and  (c)  shear 
strain  e12  illustrate  fairly  high  concentration  of 
strain  near  the  gel-kidney  interface. 
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the  product  of  <£d  and  local  displacement.  There  is  a 
limit,  however,  to  the  advantages  gained  by  increasing 
phase  shift.  As  spatial  phase  gradients  become  large,  the 
phase  distribution  within  a  given  voxel  reduces  signal 
amplitude.  Assuming  a  linear  phase  distribution  of  range 
j3  with  a  voxel,  the  signal  modulation  function  is  given  by 
sinc(j3).  For  illustration,  consider  the  vertical  phase  ex¬ 
cursion  of  21  7r  observed  across  the  80-mm  phantom  in 
Fig.  2b.  If  the  vertical  phase  excursion  was  evenly  dis¬ 
tributed  across  voxels,  the  phase  range  within  each 
0.78-mm  voxel  would  be  approximately  0.2  tt.  This  im¬ 
plies  a  signal  reduction  factor  of  sinc(0.2  tt)  =  0.94  (i.e., 
signal  loss  of  6%).  Clearly  phase  gradients  can  be  more 
concentrated  depending  on  object  geometry,  elastic  het¬ 
erogeneity,  and  deformation  geometry.  This  concentra¬ 
tion  can  lead  to  regions  of  significant  signal  loss  in  a 
manner  analogous  to  flow  dephasing  in  conventional 
MRI.  Under  such  conditions,  a  smaller  voxel  size  can 
(paradoxically)  yield  higher  signal.  In  addition,  inspec¬ 
tion  for  significant  signal  loss  within  the  displacement- 
sensitive  magnitude  image  can  identify  high-strain  re¬ 
gions  near  soft/hard  interfaces  of  a  lesion. 

Water  diffusion  in  the  presence  of  displacement  en¬ 
coding  gradients  is  another  source  of  signal  attenuation. 
We  estimate  the  signal  reduction  factor  for  freely  diffus¬ 
ing  water  was  ^0.24  in  these  experiments  (i.e.,  76% 
signal  lost  to  diffusion  effects).  Diffusion  effects  are  less¬ 
ened  by  reducing  <bd,  or  alternatively  by  shortening  TM 
without  affecting  <bd.  In  either  case,  diffusion  effects  are 
assumed  independent  of  the  deformation  state,  and 


therefore  are  ignored  in  the 
elasticity  reconstruction.  Be¬ 
cause  displacement  phase 
shift  is  the  product  of  local 
displacement  and  <&d>  the  se¬ 
lection  of  <bd  is  somewhat  ar¬ 
bitrary  as  long  as  the  applied 
differential  deformation  is  ad¬ 
equate  for  elasticity  recon¬ 
struction.  In  these  preliminary 
experiments,  the  differential 
deformation  was  <1.5  mm 
across  the  imaged  object.  For 
multi-step  acquisitions,  as 
done  here,  good  reproducibil¬ 
ity  of  deformation  is  essential. 
Significant  variation  in  defor¬ 
mation  magnitude  over  the  ac¬ 
quisition  will  lead  to  phase  in¬ 
stability,  motion-like  artifact 
in  base  images,  and  errors  that 
propagate  through  the  elastic¬ 
ity  reconstruction.  It  is  a  minor 
technical  challenge  to  achieve 
relatively  high  displacement 
reproducibility  in  the  defor¬ 
mation  apparatus.  Irreproduc- 
ible  motions  that  originate 
within  the  imaged  object, 
however,  can  be  problematic 
and  are  analogous  to  undes¬ 
ired  physiologic  motion  arti¬ 
facts  in  in  vivo  diffusion  MRI.  Fortunately,  unlike  diffu¬ 
sion  and  physiologic  motions,  the  targeted  motion  in 
elasticity  imaging  is  externally  driven.  As  such,  the  dis¬ 
placement  amplitude  in  response  to  an  external  differen¬ 
tial  deformation  can  be  significantly  greater  than  irrepro- 
ducible  or  asynchronous  displacement.  For  many  in  vivo 
applications  including  the  breast,  increasing  the  differential 
deformation  severalfold  relative  to  that  applied  in  these 
phantom  studies  can  reduce  motion  artifact.  Gradient  fac¬ 
tors  and  <t>d  would  be  reduced  accordingly,  which  would 
yield  the  added  benefit  of  increased  signal  otherwise  lost  to 
diffusion  effects. 

In  practice,  the  definition  of  a  closed  contour  of  constant 
Young’s  modulus  within  the  tissue  can  be  a  challenge.  A 
hybrid  procedure  can  be  used  as  detailed  elsewhere  (26) 
and  summarized  as  follows.  The  strain  images  are  first 
processed  to  highlight  boundaries  between  regions  of  dif¬ 
ferent  elastic  modulus.  This  procedure  is  based  on  the 
stress  continuity  property  of  continuous  media  such  as 
tissue  and  can  define  regions  of  very  small  modulus  varia¬ 
tions.  After  this  boundary  detection,  closed  contours  of 
small  elasticity  variations  are  defined.  The  modulus  along 
the  contours  is  considered  constant,  thereby  providing  the 
boundary  condition  for  complete  reconstruction  of  the  elas¬ 
ticity  within  the  region  of  interest  based  on  numerical  so¬ 
lution  of  Eq.  [6].  The  elasticity  distribution  reconstructed  in 
this  way  is  the  modulus  relative  to  the  modulus  along  the 
boundary.  For  breast  elasticity  imaging,  it  is  anticipated 
that  such  a  contour  can  be  defined  a  priori  within  the 
subcutaneous  fat  that  surrounds  breast  parenchyma.  As- 


Mapping  Young's  Modulus 


483 


FIG.  7.  (a)  Reconstructed  Young’s  modulus  image  within  a  rectangular  94  x  51  mm  region 
encompassing  the  kidney,  (b)  Relative  Young’s  modulus  for  the  scribed  ROIs  indicate  high  elastic 
modulus  at  the  giutaraldehyde  injection  site  in  the  upper-right  quadrant  of  the  kidney  parenchyma, 
(c)  The  optical  image  of  the  kidney  phantom  approximately  20  h  after  MRI.  Areas  marked  by  arrows 
were  noticeably  harder  as  assessed  by  sense  of  manual  touch. 


range  to  include  deep  lying  le¬ 
sions.  One  application  would 
be  measurement  of  elasticity 
in  breast  tissue  not  accessible 
to  manual  palpation.  In  situ 
studies  of  Young’s  elastic 
modulus  performed  on  sam¬ 
ples  of  breast  tissue  indicate 
that  there  is  a  large  difference 
in  elastic  modulus  between 
normal  and  pathologically 
transformed  breast  tissues. 
Others  have  analyzed  the 
Young’s  modulus  differences 
between  different  soft  tissues 
and  have  found  1-2  orders  of 
magnitude  difference  in 
Young’s  elastic  moduli  of  a  tis¬ 
sue  in  different  physiologic 
states  (1).  If  elastic  changes 
predate  formation  of  calcifica¬ 
tions,  elasticity  imaging  could 
potentially  increase  detection 
and/or  characterization  of  ma¬ 
lignant  breast  masses  and  thus 
be  an  important  addition  to  ex¬ 
isting  clinical  diagnostic  tools. 
Practical  issues  such  as  the  rel¬ 
atively  high  cost  of  MRI  may 
hinder  use  of  this  approach  as  a 
screening  test.  Nevertheless,  ad¬ 
ditional  work  to  define  the  role 
of  this  technique  as  a  primary 
diagnostic  tool  or  supplemental 
problem-solving  modality  in 
the  management  of  soft-tissue 
disease  is  well  justified. 


suming  the  Young’s  modulus  of  the  fat  boundary  is  con¬ 
stant,  die  relative  Young’s  modulus  image  of  parenchyma 
can  be  reconstructed.  Alternatively,  the  breast  can  be  sur¬ 
rounded  by  a  high-signal  cuff  of  known  elastic  modulus 
and  imaged.  An  image  of  absolute  Young’s  modulus  can 
then  be  reconstructed  if  the  boundary  contour  is  defined 
within  the  cuff  material. 

Some  artifacts  present  in  elasticity  images  (Figs.  4  and  7) 
are  due  to  violation  of  the  plane  strain  state  approximation 
in  these  experiments.  Indeed,  if  a  plane  strain  state  is  not 
present,  the  reconstruction  based  on  Eq.  [6]  will  be  in  error. 
The  elasticity  reconstruction,  however,  does  not  have  to  be 
limited  by  a  plane  strain  assumption  if  all  3D  components 
of  the  displacement  vector  are  available.  Fortunately,  such 
information  would  be  available  using  3D  displacement  en¬ 
coding  within  a  volumetric  imaging  sequence.  The  issue  of 
long  scan  time  could  be  resolved  by  incorporating  echo- 
planar  imaging  or  fast-spin-echo  segments  for  spatial  en¬ 
coding.  Correspondingly,  an  elasticity  reconstruction  based 
on  Eq.  [5]  would  be  applied  to  produce  volumetric  elastic¬ 
ity  maps. 

The  long-range  goal  of  quantitative  elasticity  imaging 
is  to  provide  remote  palpation  thus  expanding  its  limited 
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INTRODUCTION 

Elasticity  MRI  is  the  reconstruction  of  the  elastic  modulus  in 
media  using  local  displacement  measurements  of  an  object  under 
static  deformation  [1],  or  directly  from  measurements  of  shear 
wave  propagation  [2].  Often  it  is  assumed  out  of  plane  strain  is 
approxiamtely  zero.  However,  this  plane  strain  state 
approximation  does  not  hold  in  general.  Estimates  of  the  full  3D 
displacement  and  strain  fields  are  required  to  accurately 
reconstruct  elasticity.  Here  we  describe  such  a  method  based  on 
static  displacement,  stimulated-echo  imaging  [3].  As  previously 
described,  this  approach  encodes  local  displacement  in  response 
to  an  externally-applied  deformation.  A  relatively  long  STE 
mixing  time  (200-300msec)  allows  ill-defined  mechanical 
vibrations  to  dampen  to  where  equilibrium  conditions  apply.  The 
approach  uses  pulsed-field  gradients  (PFG)  to  encode 
displacement  and  is,  in  principle,  readily  extended  to  measure  the 
full  3D  displacement  field. 

METHODS 

Features  of  the  displacement-encoding  STE  technique  have  been 
described  previously  [3]  but  are  summarized  here.  The  initial 
spatial  configuration  of  an  elastic  object  is  encoded  ^by 
application  of  a  PFG  between  the  first  two  non-selective  90°  rf 
pulses.  An  externally-applied  deformation  force  transforms  the 
object  to  a  new  configuration  during  the  mixing  time,  TM,  of  the 
STE.  The  differential  surface  deformation  can  be  relatively 
subtle  (<5mm)  and  is  chosen  to  complement  the  applied  PFG 
area  which  defines  displacement  sensitivity  (~3-16it/mm). 
Mixing  time  must  be  sufficient  for  the  object  to  come  to  rest  by 
the  3rd  slice/slab-selective  90°  rf  pulse.  A  second  PFG  prior  to 
signal  readout  yields  a  phase  direcdy  proportional  to  local 
displacement  between  initial  and  final  object  configuration.  To 
achieve  3D,  a  fast-spin-echo  train  was  appended  to  the  STE  with 
an  additional  centric-ordered  z-slab  phase-encode  echo  train  (8- 
16  ETL  on  kz).  Conventional  phase-encoding  was  used  for  the 
remaining  spatial  dimension.  Acquisition  of  data  sets  with 
inteleaved  PFG  gradients  applied  along  X,  Y,  and  Z  axes: 
TR=l";  256  x  128  x  8;  4ave;  3.87t/mm  sensitivity.  Phase 
correction  was  done  via  a  sparse  reference  dataset  (32  ky  lines). 

A  tissue-mimicking  block  phantom  containing  ramped  bars  of 
harder  material  (=  6x  Young’s  modulus)  was  held  between  two 
pneumatic  deformation  plates.  These  parallel  plates  released  a 
5mm  differential  surface  deformation  during  TM=270ms. 

RESULTS 

Images  of  the  phantom  shown  in  figure  1  are  (a)  magnitude  and 
phase  (i.e.  displacement)  images  along  (b)  X,  (c)  Y  and  (d)  Z 
directions.  Corresponding  strain  images  illustrated  in  figure  2 

are  (a)  £„  ,  (b)  e„  ,  (c)  £xy,  and  (d)  (£»+£»)•  Note«  (e»+£yy) 
should  be  zero  if  the  plane-strain  condition  holds.  Apparent 
structure  in  figure  2(d)  indicates  there  is  out  of  plane  strain. 

DISCUSSION 

Extension  of  elasticity  MRI  to  3D  is  essential  to  implement  a 
general  elasticity  reconstruction.  3D  elasticity  reconstruction 
algorithms  have  been  described  [1],  but  have  not  been 
implemented  due  to  lack  of  suitable  3D  data  acquisition  schemes. 
The  method  presented  here  represents  an  initial  step  toward  that 
end. 


(c)  £xy  (d)  £*x+  £yy 


Figure  2 
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Abstract.  This  article  presents  a  method  for  measuring  three-dimensional  mechanical 
displacement  and  strain  fields  using  stimulated  echo  MRI.  Additional  gradient  pulses  encode 
internal  displacements  in  response  to  an  externally  applied  deformation.  By  limiting  the  mechanical 
transition  to  the  stimulated  echo  mixing  time,  a  more  accurate  static  displacement  measurement 
is  obtained.  A  three-dimensional  elasticity  reconstruction  within  a  region  of  interest  having  a 
uniform  shear  modulus  along  its  boundary  is  performed  by  numerically  solving  discretized  elasticity 
equilibrium  equations.  Data  acquisition,  strain  measurements  and  reconstruction  were  performed 
using  a  silicone  gel  phantom  containing  an  inclusion  of  known  elastic  properties.  A  comparison 
between  two-dimensional  and  three-dimensional  reconstructions  from  simulated  and  experimental 
displacement  data  shows  higher  accuracy  from  the  three-dimensional  reconstruction.  The  long¬ 
term  objective  of  this  work  is  to  provide  a  method  for  remotely  palpating  and  elastically  quantitating 
manually  inaccessible  tissues. 


I.  Introduction 

J. L  Motivation 

Palpation  has  long  been  used  by  physicians  as  a  means  to  detect  disease.  The  underlying  basis 
for  this  detection  is  the  presence  of  ‘hard’  tissue.  Evidence  suggests  that  Young’s  (or  shear) 
elastic  moduli  may  differ  by  orders  of  magnitude  within  soft  tissues  in  various  physiological 
states  (Sarvazyan  et  al  1995,  Skovoroda  et  al  1995b).  In  addition,  manual  self-examination 
is  the  first  diagnostic  line  of  defence  against  both  breast  (Hill  et  al  1988,  Newcomb  et  al 
1991)  and  testicular  cancers.  With  breast  cancer,  manual  detection  of  a  new  mass  often  merits 
excisional  biopsy,  even  if  uncorroborated  by  other  tests,  as  nodule  hardness  raises  suspicion 
of  malignancy  (Foster  1996).  Palpation  of  superficial  lymph  nodes  and  abdominal  organs  is 
also  routinely  performed.  Although  the  touch  of  a  skilled  physician  is  a  powerful  diagnostic 
tool,  palpation  sensitivity  is  relatively  poor  within  deep,  dense  or  heterogeneous  tissue.  Thus, 
most  manually  detected  lesions  are  either  superficial,  relatively  large  or  both. 

§  Author  to  whom  correspondence  should  be  addressed. 
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1.2.  Elasticity  imaging 

Currently,  many  scientists  are  working  on  extending  the  range  and  sensitivity  of  palpation  by 
using  various  methods  to  image  tissue  elasticity.  The  basic  method  for  creating  an  elasticity 
map  involves  two  steps.  First,  the  internal  displacements  within  tissue  under  an  applied 
mechanical  stress  are  measured.  The  (usually  externally)  applied  deformation  may  be  either 
dynamic  or  static.  Then,  from  these  data,  a  reconstruction  of  regional  variations  in  tissue 
elasticity  is  performed,  either  directly  or  after  calculating  internal  strains.  Although  both 
internal  displacements  and  strains  are  related  to  the  elastic  properties  of  tissue,  they  are  also 
strongly  affected  by  geometry.  Thus,  some  form  of  reconstruction  is  necessary  to  uniquely 
determine  the  elasticity  distribution. 

To  date,  two  major  medical  imaging  modalities  have  been  used  to  measure  tissue 
displacement:  ultrasound  and  magnetic  resonance  imaging  (MRI).  The  phase  sensitivity  of 
these  methods  lends  itself  to  tracking  tissue  motion.  Most  elasticity  imaging  has  been  carried 
out  using  ultrasonically  measured  tissue  displacements.  These  data  have  been  obtained  by 
tracking  specular  reflections  (Dickinson  and  Hill  1982,  Tristam  etal  1986,  1988),  by  Doppler 
techniques  (Lerner  et  al  1 990,  Parker  et  al  1990,  Parker  and  Lemer  1992),  by  cross-correlation 
of  acoustic  echoes  (Ophir  et  al  1991,  Garra  et  al  1997)  and  by  speckle  tracking  (Adler  et  al 
1989,  O’Donnell  etal  1994,  Emelianov  etal  1995).  Other  efforts  employ  MRI  for  measuring 
tissue  motion,  as  discussed  below. 


1.3.  MRI  measurement  of  tissue  displacement 

In  the  past,  myocardial  motion  and  strain  have  been  measured  using  spatial  magnetization 
tagging  (Axel  and  Dougherty  1989,  Zerhouni  etal  1988),  and  phase-based  velocity  encoding 
(Pelc  etal  1995).  More  recently,  methods  have  been  devised  to  measure  tissue  displacement 
specifically  for  elasticity  imaging.  These  measurements  can  be  separated  based  upon  the  nature 
of  the  applied  deformation. 


1.3.1.  Dynamic  deformation.  With  these  methods,  a  periodic  excitation  is  applied  to  the 
tissue  near  the  region  of  interest,  and  the  entire  system  may  be  allowed  to  reach  steady  state. 
One  or  several  ‘snapshots’  of  mechanical  wave  propagation  within  the  object  are  produced 
by  controlling  the  relative  phase  between  the  mechanical  excitation  and  the  motion-encoding 
gradients.  The  local  displacement  information  in  these  images  is  then  used  as  an  input  for 
an  elasticity  reconstruction  algorithm.  Initial  experiments  used  a  shear  excitation,  and  the 
elasticity  reconstruction  was  performed  assuming  the  recorded  image  contained  only  shear 
waves  (Muthupillai  et  al  1995).  If  only  shear  waves  are  present  in  a  purely  elastic  medium, 
local  elastic  modulus  variations  are  determined  via  the  relation  fi  =  v2X2p ,  where  p,  is 
the  local  shear  modulus,  v  is  the  frequency  of  the  applied  deformation,  X  is  the  measured 
local  strain- wave  wavelength  and  p  is  the  density  of  the  medium.  Although  attractive  in 
its  simplicity,  this  approach  is  compromised  by  frequency-dependent  viscoelastic  effects  and 
strain- wave  wavelength,  interference  from  reflections  off  of  elastic  inhomogeneities  and  the 
possible  presence  of  longitudinal  mechanical  waves  in  the  medium.  Despite  these  limitations, 
this  method  has  been  applied  in  vivo  (Dresner  etal  1999,  Lawrence  etal  1999).  Recently,  a 
more  general  elasticity  reconstruction  from  a  series  of  ‘instantaneous’  steady  state  mechanical 
wave  images  has  been  developed  (Sinkus  et  al  1999,  2000).  This  and  another  technique 
(Van  Houten  et  al  1999,  Weaver  et  al  1999)  rely  on  a  more  complete  viscoelastic  tissue  model 
than  that  presented  in  Muthupillai  etal  (1995). 
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1.3.2.  Static  deformation.  Another  method  of  producing  an  internal  strain  field  in  an  object  is 
to  deform  it  and  allow  the  material  to  relax  to  equilibrium  before  measuring  the  displacement 
field.  The  displacement  field  has  been  accessed  using  spatial  magnetization  tagging,  but  this 
method  suffers  from  spatial  resolution  limited  by  the  tagged  grid  size  and  typically  measures 
only  two-dimensional  (2D)  motion  (Fowlkes  et  al  1995).  A  quasistatic  method  using  bipolar 
gradient  phase  encoding  of  2D  motion  is  presented  by  Plewes  et  al  (1995,  1996).  Stimulated 
echo  MRI  has  also  been  used  to  measure  2D  displacement  fields  (Reese  et  al  1996),  from 
which  elasticity  images  have  been  reconstructed  (Chenevert  et  al  1998).  This  method  has 
been  extended  to  study  myocardial  motion  (Aletras  et  al  1999b).  With  these  techniques, 
viscoelastic  effects  are  generally  ignored,  making  the  reconstruction  more  straightforward. 
Care  must  be  taken,  however,  to  justify  the  use  of  a  static  model,  especially  when  repeated 
deformations  are  needed  to  acquire  a  complete  data  set. 

In  general,  MRI  has  several  advantages  over  ultrasound  with  respect  to  elasticity  imaging. 
Although  ultrasound  accurately  measures  motion  along  the  beam  axis,  lateral  motion  is 
measured  with  a  resolution  given  by  the  depth-dependent  beam  width.  Out-of-plane  motion  is 
generally  not  considered,  given  the  problems  with  three-dimensional  (3D)  image  registration 
in  ultrasound.  These  restrictions  compromise  the  quality  of  displacement  data  available  and 
constrain  the  type  of  model  used  to  produce  an  elasticity  image.  Ultrasound  does,  though, 
offer  the  advantages  of  low-cost  and  real-time  imaging.  MRI,  on  the  other  hand,  gives  one 
the  ability  to  measure  3D  displacements  within  an  object,  and  does  this  at  a  higher  overall 
resolution  than  clinical  ultrasound. 

In  this  paper  we  present  a  method  for  encoding  the  full  3D  displacement  field  within 
an  object  undergoing  an  externally  applied  static  (or  quasistatic)  deformation.  Local  strain 
estimates  are  calculated  from  the  measured  displacements,  and  the  strain  tensor  is  used  to 
numerically  solve  differential  elasticity  equilibrium  equations,  ultimately  producing  a  3D 
elasticity  image. 

2.  Reconstructive  elasticity  imaging  from  static  displacement  fields 

The  goal  of  elasticity  imaging  is  to  produce  a  map  of  the  tissue  elastic  modulus  in  a  region 
of  interest  using  available  measurements  of  displacement  components.  In  this  work,  the 
reconstruction  approach  taken  is  based  upon  a  model  of  linear,  elastic,  isotropic  media 
(Skovoroda  et  al  1995a,  1999).  The  central  equations  and  concepts  are  covered  briefly  here. 
A  more  detailed  discussion  can  be  found  in  the  references  mentioned.  Note  that  some  tissues, 
such  as  skeletal  muscle,  exhibit  anisotropic  elasticity  (Fung  1993).  For  anisotropic  media,  a 
more  generalized  reconstruction  method  is  needed. 

2. 1.  Linear  elasticity  and  reconstruction 

In  linear  elasticity,  the  components  of  the  strain  (ey)  and  stress  (oyy)  tensors  in  a  medium 
undergoing  small  deformations  are  given  by 


ffy  =  PSU  +  2Psu  (2) 


where  ut>  is  a  component  of  the  displacement  vector  U  —  (wj,  «2,  uf)  in  Cartesian  coordinates 
r  =  (x) ,  *2,  *3),  p  is  the  product  XV  •  U  for  compressible  media  or  the  static  internal  pressure 
for  incompressible  media,  Sp  is  the  Kronecker  delta  function,  X  and  p,  are  the  Lame  coefficients 
and  fi  =  pb{r)  is  the  shear  elastic  modulus. 
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A  medium  undergoing  static  deformation  obeys  the  equilibrium  condition: 


where  f  is  the  body  force  per  unit  volume  acting  in  the  jq-  direction.  In  addition,  if  a  medium 
is  incompressible,  volume  conservation  leads  to  the  following  relation: 

du\  duo 

V  •  U  =  Su  +622 +«33  =  T—  +  T—  +  T—  =0.  (4) 

0x2  ox  3 

Although  not  necessary  in  the  development  that  follows,  soft  tissue  is  approximately 
incompressible  (Sarvazyan  etal  1995). 

Using  equations  (1)  and  (2)  in  (3),  the  unknown  p(r)  can  be  eliminated  to  yield  a  set 
of  differential  equations  depending  only  on  17,  first-  and  higher-order  spatial  derivatives  of 
U,  and  the  elasticity  distribution,  This  set  of  equations  is  then  numerically  solved  to 
estimate  the  unknown  shear  elasticity  distribution. 


2.2.  Importance  of  three-dimensional  reconstruction  methods 

Several  approaches  have  been  proposed  to  estimate  tissue  elasticity  from  the  experimentally 
measured  spatial  distribution  of  internal  displacements  within  an  object.  The  simplest  method 
is  a  one-dimensional  (ID)  estimation  of  normalized  tissue  elasticity,  expressed  as 

k\  =  1/e  (5) 

where  e  is  longitudinal  strain  (Ophir  et  al  1991,  Garra  et  al  1997).  Indeed,  a  loaded  object 
generally  exhibits  low  longitudinal  strain  in  relatively  hard  regions  and  high  longitudinal  strain 
in  relatively  soft  regions. 

A  2D  elasticity  reconstruction,  based  on  a  plane-strain  assumption  and  all  necessary  in¬ 
plane  strain  components,  provides  a  more  accurate  representation  of  the  object’s  elasticity 
(Skovoroda  et  al  1995a,  1999).  The  theory  of  reconstructing  clearly  bounded  and  spatially 
distributed  tissue  inhomogeneities  has  been  demonstrated  by  Skovoroda  etal  (1995a)  as  well. 
However,  inaccurate  estimates  may  result  by  using  either  a  ID  or  2D  reconstruction  of  a  3D 
object. 

To  demonstrate  these  inaccuracies,  consider  a  spherical  inclusion  of  radius  R  in  a 
uniaxially,  uniformly  loaded,  infinite,  homogeneous  medium  (Goudier  1933).  For  an 
incompressible  medium,  the  distribution  of  longitudinal  strain  along  the  *3  axis  (orthogonal 
to  the  applied  deformation),  is  (Skovoroda  et  al  1994): 

<s 

Here  is  the  magnitude  of  the  applied  strain  and  kq  =  fi/fio  is  the  ratio  of  the  inclusion  to 
background  shear  moduli.  Normalizing  (6)  by  fi,  which  corresponds  to  the  axial  strain  in  the 
tissue  far  from  the  inclusion,  and  substituting  into  (5)  we  obtain 

3  +  2*o  R 

«>*• 

Note  that  k\/kq  =  (3  +  2ko)/5ko  within  the  inclusion.  That  is,  for  a  very  hard  inclusion 
(kq  large),  the  relative  modulus  obtained  from  a  ID  reconstruction  will  only  be  40%  of 


JC3  >  R. 


(3  +  2 Ko)  1 3  4*  2kq  4- 


*3  <  R 


x$  >  R. 
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(a)K(xrx2)  (b)K2(x1,x2) 


(e)  K3(xrx2)  (d)  Central  profiles 
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Figure  1.  Simulated  elasticity  distributions,  tc( X],xz),  and  corresponding  2D,  *2(*i.*2)  and  3D, 
K3(*i>  *2),  elasticity  reconstructions  from  the  *3 /R  ^  1.05,  ( a)-(d ),  and  *3 /R  ^  0.95,  (e)-{h), 
planes  of  a  phantom  with  a  single  hard,  spherical  inclusion  of  radius  R.  Also  presented  are  the 

central  vertical  profiles  of  each  distribution,  where  ( - )  is  k,  ( . )  is  *2 *  and  ( - )  is  *3. 

All  are  presented  on  a  log  scale  where  black  corresponds  to  a  relative  shear  modulus  of  0.5  and 
white  to  4.5.  The  background  has  a  relative  shear  modulus  of  1,  and  the  inclusion,  4. 


its  actual  value.  On  the  other  hand,  for  a  soft  inclusion  (kq  small),  the  relative  modulus 
estimate  will  approach  |,  no  matter  how  much  softer  the  inclusion  is  than  the  background. 
Obviously,  the  inaccuracy  of  a  ID  elasticity  estimation  may  not  be  acceptable  for  many 
applications. 

Now  consider  a  2D  reconstruction.  Figures  1(a)  and  ( e )  show  the  exact  relative  elasticity 
distribution,  k(x  1,  *2),  for  two  infinitesimal  planes  in  our  pedagogic  phantom  with  kq  =  4. 
Figure  1(a)  presents  k  for  x$ /R  &  1.05,  that  is,  outside  of  the  inclusion,  while  figure  \(e)  is 
the  *3 /R  &  0.95  plane.  The  corresponding  relative  2D  reconstructions,  K2(x\ ,  xf),  are  shown 
in  figures  1  (b)  and  (f ).  The  reconstructions  were  performed  using  the  algorithm  presented  by 
Skovoroda  etal  (1999).  For  comparison  with  experimental  results  (see  section  5),  an  analytic 
model  was  used  to  generate  displacement  data  which  were  sampled  with  the  *2  resolution  of  the 
experimental  displacement  encoded  data  discussed  in  section  4.2.  The  strains  used  as  input  for 
the  reconstructions  were  calculated  as  described  in  section  4.3,  and  the  reconstructions  were 
performed  over  a  region  of  interest  identical  to  the  one  discussed  in  that  same  section.  The 
positions  of  the  two  reconstructed  planes  were  selected  to  approximately  correspond  to  the 
experimental  planes  considered  in  section  5.  As  evidenced  here,  neglecting  out-of-plane  strain 
components  in  the  reconstruction  produces  geometrical  distortions  in  the  elasticity  image. 
Specifically,  the  spherical  inclusion  is  reconstructed  as  a  prolate  spheroid.  The  inaccuracy 
of  a  plane-strain  based  reconstruction  is  small  near  the  central  plane,  and  increases  with  the 
distance  between  the  imaging  plane  and  the  centre  of  the  inclusion.  Far  from  the  inclusion,  a 
2D  reconstruction  would  again  be  accurate. 
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It  is  clear  that  a  ID  or  2D  reconstruction  may  lead  to  significant  inaccuracies  in  tissue 
elasticity  estimations,  especially  when  complicated  in  vivo  geometries  influence  displacement 
and  strain  measurements.  This  points  to  the  need  for  an  accurate  3D  elasticity  reconstruction. 
A  general  unknown  shear  elasticity  distribution,  ti(x\,  x2,  *3),  must  satisfy  the  equation 
(Skovoroda  et  al  1995a,  1999): 

32(M£i2>  _  92(M£i2)  +  92[M(g22-£n)]  +  dHjie 23)  _  d2(M£i3)  _  Q 
8x I  8x\  3jCi3jC2  dx\dxi  8x28x3 

Thus,  in  order  to  compute  all  the  necessary  components  of  the  strain  tensor,  Bjj ,  in  (8),  all  of  the 
displacement  components  («i,  u2,  M3)  must  be  measured  as  a  function  of  spatial  coordinates 
(jc  1 ,  JC2 ,  *3).  This  requirement  exists  in  both  the  differential-based  3D  reconstruction  (8),  as 
well  as  in  the  more  stable  integral  based  3D  approach  (Skovoroda  et  al  1999). 

The  3D  elasticity  reconstructions  from  the  two  planes  previously  discussed  are  shown  in 
figures  1(c)  and  (g).  The  reconstruction  was  performed  as  discussed  in  section  4.3.  Although 
not  perfect  due  to  the  relatively  large  x2  step  size,  the  3D  reconstructions  clearly  exhibit 
fewer  geometric  distortions  than  the  2D  estimates.  This  is  particularly  well  illustrated  by 
the  central  vertical  profiles  through  the  analytic,  2D,  and  3D  shear  distributions  presented 
in  figures  1(d)  and  ( h ).  In  the  x2/R  ^  1.05  plane,  the  2D  reconstruction  estimates  that  an 
inclusion  is  present,  when  indeed  it  is  not,  while  the  3D  reconstruction  shows  little  evidence 
of  the  presence  of  an  inclusion.  The  estimate  of  the  extent  of  the  inclusion  in  the  x^/R  &  0.95 
plane  is  also  improved  over  the  2D  estimate.  As  with  the  2D  reconstructions,  the  strain  data 
and  reconstruction  parameters  used  for  the  3D  reconstruction  were  identical  to  those  of  the 
experimental  parameters  described  in  sections  4.2, 4.3  and  5. 

3.  Static  displacement  measurement  via  stimulated  echo  MRI 

Static  displacement  measurements  for  elasticity  imaging  avoid  several  confounding  factors  that 
may  be  present  if  dynamic  displacement  measurements  are  used.  Since  shear  wave  propagation 
speed  in  soft  tissue  is  approximately  1-20  m  s~l,  shear  waves  launched  into  a  medium  by  an 
applied  deformation  may  require  tens  of  milliseconds  to  traverse  an  object  approximately 
100  mm  in  size.  Reflected  waves  may  take  much  longer  to  dampen.  To  appropriately  measure 
an  object’s  internal  static  displacements,  the  object  must  be  in  mechanical  equilibrium — that 
is,  it  must  satisfy  (3) — during  both  the  pre-  and  post-deformation  measurements.  A  stimulated 
echo  MRI  sequence  using  displacement  encoding  gradient  pulses  is  employed  to  achieve 
this  (Reese  et  al  1996,  Chenevert  et  al  1998).  Figure  2  presents  a  schematic  of  this  pulse 
sequence.  The  mechanical  transition  from  the  pre-  to  post-deformational  states  occurs  during 
the  stimulated  echo  mixing  time,  TM.  Because  the  relevant  magnetization  is  longitudinal 
during  Tm,  it  is  unaffected  by  the  object’s  motion  during  the  mechanical  transition  period. 
This  allows  a  more  accurate  measurement  of  static  internal  displacement.  Additionally,  precise 
synchronization  of  the  motion  and  applied  gradients  is  not  necessary  as  long  as  the  mechanical 
deformation  begins  after  the  second  radio-frequency  pulse,  and  internal  motion  stops  before  the 
third.  A  long  delay  in  the  echo  time,  7e,  could  also  be  used  to  let  the  object  reach  equilibrium, 
but  this  would  likely  lead  to  prohibitive  signal  loss  from  T2  decay. 

Local  displacements  are  encoded  in  the  magnetization’s  phase  via  pulsed-field  gradients. 
The  displacement  sensitivity,  in  radians/distance,  of  the  sequence  is 

<f>d  =  y  f  Gd(r) dr  =  yGd r  (9) 

Jo 

where  y  is  the  gyromagnetic  ratio  of  the  proton,  Gd(t)  is  the  encoding  gradient  waveform  and 
r  is  the  duration  of  the  encoding  gradient.  However,  for  accurate  displacement  measurements, 
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Figure  2.  Displacement  encoding,  stimulated  echo  pulse  sequence  waveforms.  RF  =  radio 
frequency,  Gd  =  displacement  encoding  gradient,  and  Gm  =  read-out  (jcj  ),  Gpc  =  phase-encode 
fo)  and  Gsi  =  slice  (*3)  directed  gradient  waveforms.  7m  is  the  mixing  time,  7e  is  the  echo  time 
and  r  is  the  duration  of  the  displacement  encoding  gradient.  Note  that  the  displacement  encoding 
gradient  may  be  applied  to  any  of  the  directional  waveforms. 


phaseshifts  unrelated  to  the  applied  deformation  must  be  removed.  This  is  done  by  acquiring 
a  phase  reference  data  set  using  the  same  pulse  sequence,  but  with  the  object  maintained  in 
the  post-deformational  state  for  the  entire  experiment.  Note  that  all  spatial  encoding  takes 
place  with  the  object  in  the  post-deformational  state  for  both  the  displacement  encoded  and 
reference  acquisitions.  Therefore  no  image  registration  or  tracking  algorithms  are  required  to 
use  the  reference  data,  Sr,  to  correct  the  displacement  encoded  data,  Sd-  The  corrected  data 
set,  Sc,  is  then: 

£(r)  =  5d|r)5r([~r  «  |Sj(r)| e*fr>.  (10) 

\Sr(r)\ 

Most  sources  of  phase  error,  such  as  static  field  inhomogeneities,  tend  to  be  slowly  varying 
functions  of  position.  Thus  the  phase  reference  data  may  be  acquired  at  relatively  low  spatial 
resolution  to  reduce  scan  time. 

The  unwrapped  phase  of  (10)  is  related  to  the  local  displacement  vector,  ZJ,  via 

<p(r)  =  •  Ar  =  $()  •  U(r)  (11) 

where  A r  is  the  local  displacement  from  pre-  to  post-deformational  states.  The  displacement 
sensitivity,  4>d,  may  be  made  sensitive  to  motion  in  an  arbitrary  direction  based  upon 
appropriate  combination  of  displacement  encoding  gradients  in  the  read-out,  phase-encode  and 
slice  directions.  Hence,  this  pulse  sequence  readily  extends  to  acquiring  three-dimensional 
displacement  data. 

4.  Methods 

4.L  Phantom 

Elasticity  imaging  experiments  were  performed  on  a  phantom  with  a  spherical  hard  inclusion. 
Semicosil  921  silicone  gel  (Wacker  Silicones  Corporation,  Adrian,  MI)  was  used  to  construct 
a  phantom  qualitatively  simulating  the  mechanical  properties  of  soft  tissue.  The  Semicosil  921 
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consists  of  two  components,  A  and  B,  wherein  different  ratios  of  these  components  are  used 
to  vary  the  mechanical  properties  of  the  gel.  A  tissue-mimicking  phantom  was  constructed 
in  several  steps.  First,  background  material  was  prepared  by  thoroughly  mixing  components 
A  and  B  in  a  1:1  ratio,  and  then  pouring  the  mixture  into  a  154  mm  by  80  mm  rectangular 
mould.  The  mixture  was  degassed  and  cured  for  24  h  at  room  temperature  to  produce  a  22  mm 
thick  layer.  Then  a  25  mm  diameter  hard  sphere  was  prepared  from  a  1:2.5  mixture  of  A  and 
B  and  was  placed  on  top  of  the  layer  in  the  middle  of  the  mould.  Finally,  another  batch  of 
background  material  (1:1  ratio)  was  poured  into  the  mould  resulting  in  a  64  mm  by  80  mm  by 
154  mm  phantom  with  a  single,  hard,  spherical  inclusion  roughly  in  the  centre.  At  the  same 
time,  three  samples  of  each  batch  were  taken  to  independently  assess  the  elasticity  contrast 
between  the  inclusion  and  surrounding  materials.  These  measurements  were  performed  using 
the  force-deformation  system  described  in  Erkamp  etal  (1998),  and  showed  that  the  inclusion 
was  four  times  harder  than  the  background,  and  that  both  background  materials  were  elastically 
equivalent. 

4.2 .  Data  acquisition 

To  provide  repeatable  deformation,  the  phantom  was  placed  under  moderate  pre-load  pressure 
between  two  acrylic  plates  in  a  pneumatically  driven  device.  Air-filled  neoprene  boots  in  a 
push-push  configuration  provided  the  necessary  force  to  the  top  plate  to  keep  the  phantom  in 
this  pre-load  state,  and  aided  the  vertical  recoil  of  the  phantom  to  the  post-deformation  state. 
Pneumatic  pressure  was  delivered  via  two  solenoid  valves  whose  timing  was  controlled  by 
an  external  transistor-transistor  logic  circuit  triggered  by  the  pulse  sequence.  Quick-release 
valves  aided  in  depressurizing  the  boots.  Both  the  pre-load  and  recoil  positions  of  the  top 
acrylic  plate  were  set  by  adjustable  stops;  the  bottom  plate’s  position  was  fixed.  The  applied 
vertical  deformation  was  approximately  2.4  mm,  or  about  6%  strain,  between  the  pre-transition 
(greater  deformation)  and  post-transition  (less  deformation)  states. 

During  data  acquisition,  the  displacement  encoding  gradient  pulse  duration,  r ,  was  1 .5  ms, 
and  the  amplitude,  G d,  was  40  mT  m"1  in  the  read-out  (jcO  and  phase-encode  (*2)  directions, 
and  60  mT  m_1  in  the  slice  (*3)  direction.  Here,  the  x3  direction  was  along  the  bore’s  axis, 
and  the  x\  and  x2  directions  were  perpendicular  to  jc3  in  the  horizontal  and  vertical  directions 
respectively.  By  (9),  the  displacement  sensitivity,  <£d,  was  approximately  5.1 1  it  mm"1  in  the 
x\  and  jc2  directions,  and  about  7.66  n  mm-1  in  the  x3  direction.  The  displacement  encoding 
direction  was  cycled  each  pulse  repetition  between  the  jcj,  x2  and  jc3  directions.  The  pulse- 
to-pulse  repetition  time  was  approximately  0.98  s,  the  mixing  time  (FM)  was  270  ms,  and  the 
echo  time  (7E)  was  45  ms.  Two  averages  were  taken  of  a  256  x  256  x  32  matrix  covering  an 
80  mm  by  110  mm  by  48  mm  field  of  view.  The  phase  reference  data  were  collected  using 
a  256  x  32  x  32  matrix  while  keeping  all  other  parameters  the  same.  All  experiments  were 
performed  on  a  2  T,  18  cm  bore  MRI  system  (Bruker,  formerly  GE  NMR  Instruments)  using 
a  150  mm  transmit/receive  birdcage  coil. 

4.3.  Data  processing 

All  time-domain  data  were  transferred  off-line  for  processing.  For  phase  correction,  the  phase 
reference  data  set  was  zero- filled  to  a  256  x  256  x  32  matrix.  Then  this  and  the  displacement 
encoded  data  were  3D  Fourier  transformed  and  corrected  as  in  (10).  The  resulting  phase  maps 
were  then  used  to  estimate  the  spatial  derivatives  to  compute  the  strains,  via  (1),  necessary 
for  the  elasticity  reconstruction.  Phase  unwrapping  of  the  displacement  data  was  not  strictly 
required  since  only  phase  derivatives  were  used  in  the  strain  calculations.  The  displacement 
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derivative  at  the  ith  point  in  the  j  direction  was  computed  from  the  angle  of  the  complex 
multiplication  of  the  /  +  1th  point  with  the  conjugate  of  the  i  —  1th  point,  then  scaling  by 
1/20^,  where  $>Jd  is  the  magnitude  of  the  displacement  sensitivity  in  the  j  direction.  For 
convenience,  the  strain  data  were  decimated  to  the  x2  step  size  in  each  jc3  plane  in  order  to 
have  equal  resolution  in  both  the  x\  and  x2  directions.  The  strain  images  were  then  median 
filtered  with  a  5  x  5  window,  resulting  in  a  slight  decrease  in  spatial  resolution.  These  strains 
were  used  as  input  for  the  elasticity  reconstruction. 

The  3D  elasticity  reconstruction  was  performed  using  the  least-squares  error  minimization 
algorithm  discussed  in  Skovoroda  etal  (1999),  with  a  second-order,  one-sided  finite  derivative 
approximation  in  the  *3  direction.  The  reconstruction  of  fx(r)  is  a  boundary  value  problem, 
therefore  a  unique  solution  is  obtained  only  with  boundary  conditions.  So  a  square  35  mm  by 
35  mm  region  of  interest,  which  contained  the  inclusion  in  several  x3  planes,  was  identified 
in  the  x\  and  x2  directions.  Along  the  boundaries  of  these  regions,  and  in  the  two  x3  planes 
furthest  from  the  centre  of  the  inclusion  (which  did  not  contain  the  inclusion),  the  value  of  the 
shear  modulus  was  set  to  1,  resulting  in  a  relative  shear  modulus  reconstruction. 


5.  Results 

Representative  magnitude  and  corrected  phase  images  of  the  Semicosil  phantom  for  a  1 .5  mm 
thick  plane  centred  about  *3  =  0.75  mm,  or  X3/R  w  0.05,  are  shown  in  figure  3.  Knowing  that 
^5.11  7 r  mm-1  in  the  x\  and  x2  directions,  the  number  of  In  phase  wraps  in  figure  3(c) 


(a)  S1(x1,x2) 


(c)  (|)2(x1,x2) 


(b)  4>1(x1,x2) 


(d)  <t»3(xrx2) 


Figure  3.  Representative  magnitude  and  phase  images  from  the  *3 /R  &  0.05  plane  of  the  3D 
displacement  encoded  data  set  from  a  phantom  with  a  single,  hard,  spherical  inclusion.  Si,  (a),  is 
the  magnitude  of  the  *1  -displacement  encoded  data,  and  fa ,  ( b)-(d ),  are  the  phase  images  of  the 
jc,- -displacement  encoded  data. 
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Figure  4.  Representative  strain  images  from  the  X3/R  «  1.05,  (a)-(d),  and  x^/R  ^  0.95,  ( e)-Q\ ), 
planes  from  the  3D  displacement  encoded  data  of  the  phantom  with  a  single,  hard,  spherical 
inclusion.  One  normal  strain,  £22,  the  in-plane  shear  strain,  £12  =  £21,  one  through-plane  shear 
strain,fi3  =  £31 ,  and  the  trace  of  the  strain  tensor,  e\\  +£22  +  ^33,  are  presented  for  each  plane.  The 
lack  of  features  in  ( d )  and  ( h )  indicate  that  the  phantom  is  nearly  incompressible.  Linear  scales  for 
each  image  are,  from  black  to  white:  (a),  ( e ):  [-6%,  0%];  ( b ),  if ):  [-2.5%,  2.5%];  (c),  ( d ),  (#), 
(h):  [—1.6%,  1.6%]. 


indicates  a  vertical  deformation  of  approximately  2.3  mm,  and  those  in  figure  3(6)  a  horizontal 
deformation  of  about  2.0  mm.  Reduced  phase  slope  in  the  region  of  the  hard  inclusion  is 
clearly  visible  in  these  figures  as  well.  Due  to  the  central  location  of  this  plane,  there  is  little 
feature  in  <£3  (part  ( d )). 

Figure  4  shows  representative  strain  maps  from  the  planes  centred  around  *3  =  15.75  mm 
and  *3  =  14.25  mm.  Due  to  the  loaded  state  of  the  phantom  during  imaging,  the  sphere  became 
prolate,  therefore  these  planes  correspond  to  the  X3/R  &  1.05  and  X3/R  ^  0.95  locations 
respectively.  One  normal  strain,  £22  (parts  (a)  and  ( e )),  the  in-plane  shear  strain,  £12  =  £21 
(parts  ( 6 )  and  (/*)),  and  one  through-plane  shear  strain,  £13  =  £31  (parts  (c)  and  ( g )),  are  shown 
for  each  plane.  These  components  are  all  required  to  perform  the  elasticity  reconstruction  in 
(8).  Note  that  the  presence  of  through-plane  strains  in  (8)  necessitates  measurement  of  the  full 
3D  displacement  field.  In  addition,  although  elasticity-specific  details  are  seen  in  the  strain 
maps,  features  related  to  geometry  and  the  applied  deformation  are  also  clearly  present.  This 
points  to  the  need  for  a  proper  elasticity  reconstruction  to  disentangle  these  factors.  Also 
shown  is  the  trace  of  the  strain  tensor,  £11  +  £22  +  £33,  for  each  plane  (parts  (d)  and  ( h )).  The 
relative  lack  of  features  in  the  trace  of  the  strain  tensor  indicates  that  the  phantom  is  nearly 
incompressible  (like  soft  tissue). 

Magnitude  images  of  the  35  mm  by  35  mm  regions  of  interest  in  the  same  two  planes,  along 
with  two  different  shear  modulus  reconstructions  of  these  planes,  are  presented  in  figure  5. 
In  the  magnitude  images,  the  hard  inclusion  is  clearly  present  in  the  x$/R  «  0.95  plane, 
while  it  is  essentially  absent  in  the  x^/R  &  1.05  plane.  Note  that  the  magnitude  images  only 
convey  geometric  information.  Figures  5(6)  and  (f)  show  2D  elasticity  reconstructions  of  these 
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(f)  K2(xrx2) 


(g)  k3(x1?x2) 
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Figure  5.  x\ -displacement  encoded  magnitude  images,  SiUi,  *2),  and  corresponding  2D, 
*2(*i**2).  and  3D,  ^3 ( jc j ,  JC2 )  elasticity  reconstructions  from  the  *3 /R  «  1.05,  ( a)-(d ),  and 
xjj R  &  0.95,  ( e)-(h ),  planes  of  a  phantom  with  a  single  hard,  spherical  inclusion.  Also  presented 

are  the  central  vertical  profiles  of  each  distribution,  where  ( . )  is  kj,  and  (----)  is  *3.  All 

are  presented  on  a  log  scale  where  black  corresponds  to  a  relative  shear  modulus  of  0.5  and  white 
to  4.5.  The  background  has  a  relative  shear  modulus  of  1,  and  the  inclusion,  4,  from  independent 
measurements.  For  geometric  reference  purposes,  parts  ( d )  and  ( h )  include  plots  of  Sj"1  (filtered 
and  normalized)  as  ( - ). 


two  planes,  while  figures  5(c)  and  ( g )  show  the  corresponding  3D  elasticity  reconstructions. 
As  in  figure  1,  one  sees  an  overestimate  of  the  3D  spatial  extent  of  the  inclusion  in  the 
2D  reconstructions.  This  overestimate  is  corrected  with  the  3D  reconstruction.  For  ease 
of  comparison,  vertical  profiles  through  the  centre  of  the  inclusion  from  the  2D  and  3D 
reconstructions  are  presented,  along  with  plots  of  Sj-1  for  geometric  reference,  in  figures  5(d) 
and  (h). 


6.  Discussion 

The  stimulated  echo  sequence  presented  here  phase  encodes  internal  displacements  using 
gradient  pulses.  An  externally  applied  deformation,  synchronized  with  the  pulse  sequence, 
produces  an  internal  displacement  field.  This  deformation  is  actively  driven  with  a  pneumatic 
device,  and  the  mechanical  transition  from  pre-  to  post-deformation  occurs  during  the  sequence 
mixing  time,  7m,  while  the  relevant  magnetization  is  longitudinal.  Because  longitudinal 
magnetization  decays  only  as  T\,  the  mechanical  transition  period  may  be  extended  to  allow 
potentially  long-lived  or  ill-defined  motions  within  the  object  to  dampen.  With  a  sufficiently 
long  7m,  the  encoded  displacement  will  be  approximately  static.  However,  signal  loss  due  to 
T\  relaxation  sets  a  practical  limit  on  the  length  of  TM.  To  determine  an  appropriate  mixing 
time,  a  series  of  2D  displacement  encoded  images  was  taken,  varying  rM  from  50-750  ms. 
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From  the  phase  maps  of  these  data,  one  could  see  that  the  top  acrylic  plate  of  the  deformation 
device  completed  its  excursion  in  under  200  ms,  and  the  internal  motion  in  the  phantom 
became  negligible  by  250  ms.  To  ensure  that  (3)  was  reasonably  satisfied,  a  270  ms  mixing 
time  was  chosen  for  subsequent  data  collection.  Although  the  deformation  device  used  here 
provides  adequate  transition  speed,  an  even  faster  device  would  allow  a  shorter  TM,  yielding 
more  signal.  While  absent  in  the  phantom  used  here,  water  diffusion  in  the  presence  of 
displacement  encoding  gradients  will  be  another  source  of  signal  loss  in  in  vivo  experiments. 
This  loss  can  be  mitigated  by  reducing  the  displacement  sensitivity,  <J>j,  or  by  shortening  7m- 

The  quality  of  the  shear  modulus  reconstruction  ultimately  depends  on  the  local  phase, 
as  in  (1 1),  induced  by  the  encoding  gradients  and  the  local  displacement.  More  specifically, 
the  quality  depends  on  the  spatial  derivatives  of  the  encoded  phase.  A  study  of  the  effects 
of  displacement  sensitivity,  applied  deformation,  relative  hardness  and  diffusion  loss  on  the 
signal-to-noise  ratio  (SNR)  of  the  elasticity  reconstruction  has  been  presented  in  Steele  et  al 
(1999).  This  study  demonstrates  that  increased  intra- voxel  phase  wrap  will  increase  the 
reconstruction  SNR,  up  to  a  tt  intra- voxel  phase  distribution.  Note  that  the  reconstruction 
SNR  increases  despite  a  reduction  in  the  nuclear  magnetic  resonance  (NMR)  signal  from  the 
object.  Assuming  a  linear  phase  distribution  of  0  radians  across  a  voxel,  the  signal  modulation 
from  that  voxel  will  be  |sinc(0/2)|  =  |  sin(0/2)(0/2)-1|.  However,  the  phase  gradients  (that 
is,  the  displacement  derivatives)  will  be  maximized  without  aliasing  as  the  intra-voxel  phase 
wrap  approaches  tt,  and  this  is  the  signal  that  is  important  in  the  reconstruction.  A  tv  phase 
wrap  may  be  achieved  through  many  combinations  of  applied  deformation  and  displacement 
sensitivity.  However,  increasing  4>d  will  increase  signal  loss  due  to  diffusion,  as  discussed 
above.  Hence,  a  smaller  displacement  sensitivity  and  increased  deformation  would  appear  to 
be  optimal.  Again,  there  is  a  trade-off:  as  deformation  increases,  the  model  of  linear  elasticity 
discussed  in  section  2.1  will  become  less  and  less  valid.  Elasticity  reconstructions  from  finite 
displacement  fields  have  been  demonstrated  in  Skovoroda  et  al  (1999),  but  these  are  obviously 
more  computationally  intensive  than  the  linear  reconstructions  used  here.  In  relation  to  the  data 
presented  here,  the  number  of  27T  phase  bands  across  the  phantom  in  figure  3  clearly  indicate 
that  these  data  were  acquired  with  a  suboptimal  displacement  sensitivity/applied  deformation 
combination.  Because  neither  the  encoding  nor  the  deformation  used  here  were  extreme,  the 
elasticity  reconstruction’s  SNR  should  be  improved  merely  by  optimizing  the  intra-voxel  phase 
wrap. 

Relative  hardness,  object  geometry  and  deformation  geometry  also  affect  the  displacement 
phase  gradients.  In  general,  the  phase  gradients  increase  near  soft/hard  interfaces  and  are  higher 
in  relatively  soft  regions  of  tissue.  Excessive  phase  wrap  (i.e.  strain)  can  lead  to  regions  of 
significant  signal  loss  in  a  manner  analogous  to  flow  dephasing  in  conventional  MRI.  The 
resulting  reconstructions  would  suffer  from  this  signal  loss.  Hence,  the  applied  deformation 
and  displacement  sensitivity  should  be  optimized  for  the  regions  of  highest  strain  in  an  object. 
Increased  intra- voxel  phase  wrap  in  regions  of  lower  strain  may  be  obtained  by  integrating  the 
signal  from  several  voxels;  in  essence,  applying  an  adaptive  voxel  size  based  upon  local  phase 
gradients.  This  increase  in  signal  would  come  at  the  expense  of  spatial  resolution.  Additionally, 
signal  loss  due  to  intra-voxel  phase  wrap  in  the  displacement  encoded  magnitude  images  may 
be  useful  for  identifying  regions  of  high  strain  in  tissue. 

Another  factor  affecting  the  displacement  signal  is  the  reproducibility  of  the  applied 
deformation.  For  multistep  acquisitions,  such  as  those  presented  here,  good  deformation 
reproducibility  is  essential.  Variations  in  the  applied  deformation  will  lead  to  phase  instability, 
motion-like  artefacts,  and  errors  that  will  propagate  through  the  elasticity  reconstruction. 
Adequate  reproducibility  has  been  achieved  with  the  current  deformation  system.  However, 
irreproducible  or  asynchronous  motions  within  the  imaged  object  may  be  problematic. 
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These  would  include  physiological  cardiac  and  respiratory  motion  present  in  in  vivo 
experiments.  In  some  ways,  the  problems  associated  with  undesired  motion  would  be  similar 
to  those  encountered  in  diffusion  MRI.  Because  the  applied  deformation  is  external,  though, 
the  displacement  encoding  can  be  tailored  to  it,  reducing  the  effect  of  undesired  motion  on  the 
displacement  data.  Further  complications  arise  because  phase  derivatives  of  the  displacement 
data,  approximated  by  finite  differences,  are  required  for  reconstruction.  In  addition  to 
choosing  an  appropriate  deformation/encoding  combination,  methods  should  be  devised  to 
reduce  the  effects  of  undesired  motion  on  the  displacement  derivatives. 

Clearly  several  advantages  justify  performing  a  3D  elasticity  reconstruction  rather  than  a 
2D  reconstruction.  As  illustrated  in  figures  1  and  5,  and  as  discussed  in  sections  2.2  and  5,  a 
3D  reconstruction  provides  a  more  accurate  representation  of  the  elasticity  distribution  than 
a  2D  reconstruction  in  the  simple  phantom  used  here.  Complicated  in  vivo  geometries  will 
only  increase  the  likelihood  that  neglecting  out-of-plane  strain  components  will  result  in  an 
inaccurate  elasticity  estimate.  This  increased  accuracy  comes,  though,  at  the  expense  of 
increased  computational  complexity  and  increased  scan  time.  For  instance,  a  single  acquisition 
of  the  3D  data  discussed  in  section  4.2  takes  over  6  h!  This  far  exceeds  any  clinically  feasible 
scan  times.  The  total  experiment  time  may  be  lessened  through  the  use  of  echo-planar  imaging 
(Mansfield  1977)  or  fast  spin-echoes  (Hennig  etal  1986)  for  spatial  encoding  (Chenevert  etal 
1999).  A  fast  scan  implementation  has  already  been  used  to  study  cardiac  motion  (Aletras  etal 
1999a).  The  number  of  planes  of  data  acquired  may  also  be  significantly  reduced  while  still 
allowing  a  3D  reconstruction,  shortening  the  scan  time  further.  It  should  be  mentioned  that  a 
classic  3D  stimulated  echo  sequence  was  deliberately  chosen  in  part  due  to  SNR  considerations, 
since  a  fast  scan  implementation  of  the  method  would  generally  have  a  lower  SNR  than  one 
classically  phase  encoded.  Being  an  inverse  problem,  the  3D  reconstruction  is  sensitive  to  the 
SNR,  and  we  wanted  the  initial  test  of  the  reconstruction  to  be  done  with  the  highest  SNR  data 
possible  using  this  technique.  Also,  note  that  the  reconstruction  in  (8)  does  not  rely  on  the 
assumption  of  incompressibility,  although  making  this  assumption  provides  another  means  of 
regularizing  the  inverse  problem. 

Additionally,  a  reconstruction  of  static  displacement  data  offers  several  advantages  over 
a  reconstruction  of  dynamic  displacement  data.  A  static  reconstruction  allows  one  to  ignore 
viscoelastic  effects  as  well  as  the  longitudinal  or  shear  nature  of  the  applied  deformation.  Static 
methods  also  provide  high  SNR  displacement  and  strain  estimates.  Dynamic  methods,  on  the 
other  hand,  provide  a  potentially  very  simple  reconstruction  (Muthupillai  etal  1995).  However, 
this  reconstruction  may  be  compromised  by  interference  from  elastic  inhomogeneities, 
attenuation  of  shear  waves,  mixing  of  longitudinal  and  shear  waves,  and  resolution  limits 
imposed  by  noise  when  determining  the  shear- wave  wavelength.  Reconstruction  models  that 
include  viscoelastic  effects  allow  a  more  accurate  interpretation  of  dynamic  data  (Sinkus  et  al 
1999,  2000,  Van  Houten  et  al  1999),  but  these  are  necessarily  more  complicated  than  static 
models  (Skovoroda  et  al  1995a,  1999). 

Choosing  a  contour  of  constant  shear  modulus  for  appropriate  boundary  conditions  for  (8), 
though,  can  in  practice  be  a  challenge.  In  the  applications  discussed  here,  a  priori  knowledge 
of  phantom  geometry  was  employed  in  the  reconstructions.  This  may  be  possible  in  vivo  as 
well,  albeit  more  complicated.  For  instance,  in  breast  elasticity  imaging,  such  a  contour  may 
be  defined  in  the  subcutaneous  fat  surrounding  the  parenchyma  using  the  boundary  detection 
procedure  described  in  Skovoroda  et  al  (1995a).  The  elasticity  reconstruction  would  then  be 
relative  to  the  shear  modulus  of  the  fat  boundary,  assuming  that  it  is  constant.  Alternatively,  a 
high  signal  cuff  of  known  elastic  modulus  could  be  used  to  surround  the  breast.  This  would 
provide  an  absolute  image  of  shear  modulus  variations  if  the  boundary  contour  were  chosen 
inside  the  cuff. 
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The  3D  shear  elasticity  reconstructions  presented  above  contain  artefacts  both  inside 
and  outside  of  the  inclusion  due  to  the  finite  SNR  in  the  measured  displacement  strain 
components,  and  due  to  the  step  size  used  in  the  finite  approximation  to  the  derivatives  in  the 
reconstruction.  In  contrast  to  the  2D  elasticity  reconstruction,  where  the  elasticity  distribution  is 
reconstructed  independently  in  each  plane,  the  3D  reconstruction  uses  the  elasticity  distribution 
in  neighbouring  planes.  Therefore,  in  addition  to  in-plane  error  propagation  problems 
discussed  elsewhere  (Skovoroda  et  al  1995a),  error  propagation  in  the  through-plane  direction 
may  occur  due  to  inaccurate  elasticity  reconstructions  in  the  preceding  planes.  This  is 
particularly  true  if  the  3D  elasticity  reconstruction  is  performed,  as  in  this  paper,  by  solving  an 
initial  value  problem  in  the  through-plane  direction.  Even  though  the  more  stable  integral  based 
approach  (Skovoroda  et  al  1999)  was  employed  to  solve  for  fi(r)  in  each  plane,  the  results 
of  the  3D  elasticity  reconstructions  in  subsequent  planes  exhibit  significant  error  propagation 
in  the  *3  direction.  Given  a  particular  spatial  discretization  of  the  displacement  data,  this 
error  propagation  can  be  reduced  by  several  approaches.  These  include  more  appropriate  data 
filtering  and  reducing  the  reconstruction’s  sensitivity  to  noise,  but  these  considerations  are 
beyond  the  scope  of  this  paper. 

7.  Conclusions 

The  ultimate  goal  of  quantitative  elasticity  imaging  is  to  provide  physicians  with  a  method 
of  remotely  palpating  soft  tissue  to  detect  disease.  The  three-dimensional  elasticity  imaging 
technique  demonstrated  here  is  a  step  toward  extending  the  range  and  sensitivity  of  palpation, 
a  powerful  diagnostic  tool.  One  possible  application  of  this  technique  would  be  measuring 
the  elasticity  of  breast  tissue  normally  inaccessible  to  manual  palpation.  A  large  elastic 
modulus  difference  between  normal  and  pathological  breast  tissue  has  been  measured  in  situ. 
A  previous  study  indicates  that  soft  tissues  in  different  physiological  states  display  shear 
modulus  variations  of  one  to  two  orders  of  magnitude  (Sarvazyan  et  al  1995).  If  these  elastic 
changes  pre-date  calcification  formation,  elasticity  imaging  may  increase  sensitivity  to  and 
characterization  of  malignant  breast  masses,  complementing  existing  diagnostic  tools.  The 
relatively  high  cost  of  MRI  may  hinder  using  this  approach  as  a  general  screening  technique. 
However,  additional  work  to  define  the  role  of  this  modality  as  a  primary  or  complementary 
diagnostic  tool  in  diseases  of  soft  tissues  seems  worthwhile  indeed. 
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INTRODUCTION 

Elasticity  imaging  may  offer  clinicians  a  non-invasive 
method  of  remotely  palpating  patients  to  detect  diseased  tissues. 
Tissue  elasticity  information  may  be  obtained  via  NMR  detection 
of  shear  waves1  or  of  static  tissue  displacement.2  Here  we 
present  simulation  results  elucidating  the  effects  of  displacement 
sensitivity  and  applied  deformation  on  the  signal-to-noise  ratio 
(SNR)  of  elasticity  images  obtained  using  static  displacement, 
stimulated  echo  NMRI  (SDSEI). 

PRINCIPLE 

In  SDSEI,  the  object  is  deformed  during  the  stimulated  echo 
mixing  time,  TM.  Local  displacements  are  encoded  by  means  of 
a  pulsed  gradient  applied  before  and  after  deformation.  The 
sensitivity  is  given  by 

*d=rGdT 

where  the  displacement  sensitivity,  Od,  is  in  radians/distance,  Gd 
is  the  displacement  encoding  gradient  strength,  and  T  is  the 
displacement  encoding  gradient  duration.  The  gradient  of  the 
displacement  field  yields  the  strain  field,  and  the  strain  field  is 
used  as  an  input  for  a  boundary  value  problem  to  extract  the 
tissue  elastic  moduli.3 
SIMULATION 

A  simple  one-dimensional,  100-mm  long  object  is  used  in 
all  simulations.  This  object  contains  an  inclusion  whose  shear 
modulus  is  five  times  larger  than  the  surrounding  material’s. 

The  pulse  sequence  assumed  is  that  in  the  paper  by  Chenevert  et 
a/.,2  with  TM  =  200  ms,  x  =  4.5  ms,  and  225  ms  between 
displacement  encoding  gradients  .  The  NMR  resolution  is  128 
pixels.  The  base  SNR  of  the  material  is  set  to  130,  matching 
the  SNR  of  a  tissue-mimicking  material  measured  using  the 
SDSEI  pulse  sequence. 

The  variables  used  are  Gd  and  the  normalized  applied  surface 
deformation,  Go-  Since  x  is  assumed  constant,  Od  varies  linearly 
with  Gd.  Gd  ranges  from  0  to  5  G/cm  in  0.2  G/cm  increments, 
and  Go  from  1%  to  20%  in  1%  steps.  For  each  gradient- 
deformation  combination,  both  the  NMR  image  and  the  elasticity 
image  are  reconstructed.  The  reconstructed  elastic  moduli  are 
normalized  to  the  background  reconstruction.  Then  the  NMR 
image  SNR,  SNRo,  and  the  elasticity  image  SNR,  SNRE,  are 
calculated.  Four  different  noise  instances  are  averaged  for  every 
point. 

RESULTS 

Figures  la  and  lb  show  the  SNR0  and  the  SNRE, 
respectively,  of  the  background.  Here  we  see  that  the  SNRo 
decreases  with  increasing  Gd  and  Gq  due  to  intrapixel  dephasing. 
Note  that  signal  loss  due  to  diffusion  increases  as  Gd  increases. 
The  SNRe  begins  at  0  for  Gd  =  0,  and  then  increases  with  both 
deformation  and  sensitivity  until  the  intrapixel  dephasing  is  n 
radians.  Beyond  this  point  it  drops  off  quickly  to  zero.  Again, 
diffusion  signal  loss  is  present. 

Figures  2a  and  2b  show,  respectively,  the  SNRo  and  die 
SNRe  of  the  inclusion.  Here  the  major  source  of  signal  loss  in 
SNR0  is  diffusion.  Because  the  inclusion  is  harder  than  the 
background,  it  tends  to  displace  rather  than  deform  under  the 
applied  surface  deformation;  hence  only  the  gradient  contributes 
significantly  to  the  intrapixel  dephasing.  The  SNRE  simply 
seems  to  exhibit  diffusion  loss  and  does  not  peak  as  the  SNRE  of 
the  background  does.  Also  note  that  the  SNRE  of  the  inclusion 
is  lower  than  that  of  the  background. 


This  can  be  explained  by  noting  that  the  elasticity 
reconstruction  depends  upon  the  strain  field,  and  the  strain  is  the 
gradient  of  the  displacement  field.  As  noted  above,  the  hard 
inclusion  tends  to  displace,  not  deform,  under  a  certain  surface 
deformation.  Thus,  not  only  does  the  intrapixel  dephasing 
decrease,  but  the  m/erpixel  dephasing  (i.e.,  strain)  decreases. 
Hence  the  combinations  of  Od  and  g0  presented  here  do  not 
approach  the  optimal  strain  fields  for  the  inclusion 
reconstruction. 


Fig.  la:  Background  SNR0 


Fig.  2a:  Inclusion  SNRo 
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CONCLUSIONS 


We  have  illustrated  the  effect  of  displacement  sensitivity  and 
applied  deformation  on  static  displacement  NMR  elasticity 
images  and  noted  the  differing  optimal  imaging  parameters  for 
tissues  with  different  shear  moduli.  These  results  suggest  a  type 
of  adaptive  elasticity  reconstruction  wherein  the  voxel  size  for  a 
given  region  of  tissue  varies  inversely  with  the  strain  field  in 
order  to  optimize  the  elasticity  reconstruction’s  SNR. 
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Abstract  A  method  is  presented  to  reconstruct  the  elas¬ 
tic  modulus  of  soft  tissue  based  on  ultrasonic  displacement 
and  strain  images  for  comparatively  large  deformations.  If 
the  average  deformation  is  too  large  to  be  described  with 
a  linear  elastic  model,  nonlinear  displacement-strain  rela¬ 
tions  must  be  used  and  the  mechanical  equilibrium  equa¬ 
tions  must  include  high  order  spatial  derivatives  of  the  dis¬ 
placement.  Numerical  methods  were  developed  to  reduce 
error  propagation  in  reconstruction  algorithms,  including 
these  higher  order  derivatives.  Problems  arising  with  the 
methods,  as  well  as  results  using  ultrasound  measurements 
on  gel-based,  tissue  equivalent  phantoms,  are  given.  Com¬ 
parison  to  reconstructions  using  a  linear  elastic  model  shows 
that  equivalent  image  quality  can  be  produced  with  algo¬ 
rithms  appropriate  for  finite  amplitude  deformations. 


I.  Introduction 

TMages  of  mechanical  displacements  and  strains  within 
Xsoft  tissue  present  information  about  the  elasticity  of  in¬ 
ternal  structures  [1]— [21].  Interpreting  these  images,  how¬ 
ever,  can  be  difficult  for  complex  mechanical  objects  such 
as  soft  tissue.  To  potentially  simplify  image  interpretation 
and  reduce  artifacts  due  solely  to  object  geometry,  several 
investigators  have  explored  elastic  modulus  reconstruction 
[9],  [10],  [17],  [22].  Exact  reconstruction  is  impossible  with- 
out  detailed  knowledge  of  the  mechanical  boundary  condi¬ 
tions  (i.e.,  the  Young’s  modulus  needs  to  be  specified  along 
some  boundary,  as  discussed  later  in  this  paper).  Never¬ 
theless,  methods  have  been  developed  to  produce  relative 
reconstructions,  even  if  detailed  mechanical  boundary  con¬ 
ditions  are  unknown  [9],  [10]. 

Previous  work  from  our  laboratories  has  shown  that  the 
Young’s  (or  equivalently  the  shear)  modulus  of  soft  tissue 
and  tissue-like  phantoms  can  be  reconstructed  from  me¬ 
chanical  displacement  and  strain  images  acquired  during 
static  external  deformation  [5]-[10],  [12],  [18]-[21],  Differ¬ 
ent  imaging  systems  (e.g.,  ultrasound  and  MRI),  as  well 
as  different  deformation  procedures,  were  used  to  gener¬ 
ate  displacement  and  strain  images.  Consequently,  effec- 
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tive  numerical  methods  were  developed  for  all  systems  to 
reconstruct  the  relative  Young’s  modulus  based  on  a  linear 
elastic  model.  These  techniques  do  not  require  any  infor¬ 
mation  about  global  boundary  conditions  (i.e.,  mechanical 
constraint  of  the  body,  its  geometry,  the  types  of  external 
and  internal  forces,  etc.  [9]).  In  principle,  however,  they  are 
limited  to  low  magnitude  external  deformations  in  which 
a  linear  model  is  valid.  Here  we  extend  these  methods  to 
finite  amplitude  deformations. 

Large  external  deformations  increase  the  signal-to-noise 
ratio  (SNR)  of  displacement  and  strain  images  [5]-[10], 
[23],  Unfortunately,  large  deformations  of  soft  tissue  and 
tissue-like  materials  cannot  be  described  with  a  linear  elas¬ 
tic  model.  A  linear  model  can  break  down  in  two  ways. 
First,  for  most  soft  tissues,  the  elastic  modulus  increases 
as  a  function  of  strain  (i.e.,  strain  hardening).  This  effect 
is  often  referred  to  as  “material  nonlinearity.”  Second,  a 
more  complete  description  of  the  equilibrium  equation,  in¬ 
cluding  nonlinear  strain-displacement  relations,  must  be 
used  for  large  deformations.  This  effect  is  often  referred 
to  as  “geometric  nonlinearity.”  Due  to  the  high  order  dis¬ 
placement  derivatives  resulting  from  this  description,  error 
propagation  must  be  minimized  in  any  reconstruction  al¬ 
gorithm  using  measured  displacement  data  with  a  finite 
signal  to  noise  ratio. 

In  this  paper,  we  examine  the  second  form  of  nonlin¬ 
earity,  namely,  geometric  nonlinearities  arising  from  large 
amplitude  deformations.  These  studies  were  conducted  on 
gelatin  phantoms  with  almost  no  material  nonlinearities 
over  the  deformation  range  considered  here  (average  strain 
up  to  14%).  Material  nonlinearities  in  soft  tissue  are  con¬ 
sidered  in  [24],  [25].  The  specific  purpose  of  the  present 
study  was  to  explore  numerical  methods  minimizing  the 
effects  of  higher  order  displacement  derivatives  needed  to 
describe  finite  amplitude  deformations  on  elasticity  recon¬ 
struction. 

Previous  algorithms  for  elasticity  reconstruction  were 
formulated  using  the  set  of  equations  describing  mechan¬ 
ical  equilibrium  in  a  statically  deformed,  linear  elastic 
medium  [5],  [9],  [10],  [18]-[21],  [24],  [26],  Independent  of 
the  specific  elastic  model,  however,  these  equations  can  be 
posed  in  either  differential  or  integral  form.  An  integral 
representation  is  more  appropriate  for  a  nonlinear  model 
given  realistic  measurement  noise.  As  discussed  in  Sec¬ 
tion  II,  numerical  methods  have  been  developed  for  both 
linear  and  nonlinear  models  exploiting  an  integral  repre¬ 
sentation  of  the  reconstruction  equations.  The  specific  ap¬ 
proach  assumes  a  plane  strain  state  to  approximate  two- 
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dimensional  displacement  and  strain  images  obtained  with 
a  real-time  ultrasound  scanner.  The  full  three-dimensional 
problem  is  discussed  in  the  Appendix. 

Displacement  data  acquired  with  a  real-time  ultrasound 
scanner  were  used  to  test  the  numerical  methods  of  Sec¬ 
tion  II.  Relative  elastic  modulus  images  were  reconstructed 
within  a  gel-based,  tissue  equivalent  phantom  with  pre¬ 
scribed  mechanical  properties  using  both  linear  and  non¬ 
linear  models.  Methods  for  displacement  and  strain  image 
acquisition  are  presented  in  Section  III.  All  results  are  pre¬ 
sented  in  Section  IV.  The  paper  concludes  with  a  discus¬ 
sion  of  the  results  in  Section  V. 


II.  Theory 

Consider  a  three-dimensional  (3-D)  volume  V  of 
deformed  media  with  the  displacement  vector  U  = 
U(x  1,^2, X3)  =  {uuU2,Us)  in  Cartesian  coordinates  rr*, 
i  =  1,2,3.  The  volume  V  can  be  either  the  entire  mechan¬ 
ical  body  or  a  region  of  interest  inside  the  object. 

The  most  general  nonlinear  mechanical  equilibrium 
equations  are  [27]-[31]. 

E  {&»#<»  +  «<,»)])  =0  *  =  1-2,3.  (1) 

j—l  [n=l  J  j 

Here  crnj  is  a  component  of  the  2nd  ranked  stress  ten¬ 
sor  and  din  is  the  Kronecker  delta  symbol.  In  (1),  and  the 
rest  of  this  paper,  the  lower  index  after  a  comma  means 
differentiation  with  respect  to  the  corresponding  spatial 
Lagrangian  coordinate.  Note  that  spatial  coordinates  and 
displacement  components  correspond  to  the  initial,  not  de¬ 
formed,  configuration  of  the  object  under  investigation. 
Similarly,  all  images  are  presented  in  the  original  object 
geometry  (i.e.,  before  deformation). 

Equation  (1)  must  be  satisfied  at  every  internal  point 
of  the  body.  If  the  magnitudes  of  the  spatial  derivatives 
of  all  displacement  components  are  small,  the  last  terms 
UiiTl  in  (1)  can  be  omitted,  producing  the  familiar  linear 
equilibrium  equations  [32],  [33]: 

3 

E^,i  =  0,  i  =  1,2,3.  (2) 

3= 1 

To  complete  the  system  of  equations  describing  internal 
deformations,  the  relation  between  stress  and  strain  ten¬ 
sors,  as  well  as  the  relation  between  the  strain  tensor  and 
the  displacement  vector,  are  needed.  Here  we  assume  that 
the  standard  linear  relation  between  the  stress  tensor  err¬ 
and  the  strain  tensor  Eij  for  incompressible  media  is  still 
valid  [27]— [34] : 

&ij  —  pdij  T  2{i£ij ,  (3) 

where  p  is  the  static,  internal  pressure  and  the  shear  elas¬ 
tic  modulus  p  is  considered  a  constant  independent  of  the 
strain  magnitude.  Computing  the  spatial  distribution  of 


the  shear  elastic  modulus  is  the  goal  of  reconstruction. 
Note  in  an  incompressible  material,  such  as  soft  tissue, 
the  shear  and  Young’s  moduli  are  simply  proportional  (i.e., 
E  =  /i/3).  Thus,  shear  modulus  and  Young’s  modulus  re¬ 
constructions  are  equivalent. 

In  the  following  derivation  we  assume  a  plane  strain 
state  in  which  spatial  derivatives  of  the  out-of-plane  dis¬ 
placement  us  are  either  zero  or  small  compared  to  the 
others,  and  the  two  in-plane  components  U\  and  u 2  do  not 
vary  significantly  as  functions  of  the  out-of-plane  coordi¬ 
nate  [5],  [6],  [8],  [10],  [12],  [18],  [19],  [21],  [32].  Nonlinear 
elasticity  reconstruction  for  a  general  three-dimensional 
strain  state  is  considered  and  discussed  in  the  Appendix. 

Because  the  pressure  p  cannot  be  directly  measured 
with  an  imaging  system,  it  must  be  eliminated  from  the 
equations  describing  mechanical  equilibrium  [5],  [9],  [10], 
[18]— [22] ,  [26].  For  the  linear  case,  eliminating  p  from  (2) 
and  (3)  leads  to  a  partial  differential  equation  for  the  un¬ 
known  shear  elasticity  distribution  y,  —  /x(x,  y): 

(p£xy),xx  —  {p£xy),yy  H"  %{P'£yy) ,xy  =  0?  (4) 

where  the  spatial  coordinates  x\  and  X2  are  denoted  here 
by  x  and  y.  Note  that  the  incompressibility  condition,  £n  + 
622  =  0,  is  used  to  produce  the  specific  form  presented  in 
(4).  All  strains  in  this  equation  are  defined  by  the  linear 
strain-displacement  relation: 

4jn  =  \(^j  +Uj,i)-  (5) 

The  reconstruction  procedure  assumes  that  all  relevant 
displacements  are  known  (i.e.,  measured),  and  solves  for 
the  shear  elasticity  distribution  satisfying  both  (4)  and 
the  mechanical  boundary  conditions. 

Except  for  the  degenerate  (i.e.,  parabolic)  case  corre¬ 
sponding  to  an  in-plane  translation  with  rotation  of  the 
volume  as  a  rigid  body,  internal  deformations  are  described 
by  a  hyperbolic,  2nd  order  differential  equation  [9],  [26]. 
This  means  there  is  a  unique  shear  modulus  distribution 
satisfying  (4)  given  appropriate  mechanical  boundary  con¬ 
ditions.  Rather  than  specifying  displacement  and  stress 
values  at  the  object  boundary  [17],  a  much  simpler  bound¬ 
ary  value  problem  can  be  formulated  using  the  method 
of  characteristics  in  which  the  modulus,  and/or  spatial 
derivatives  of  the  modulus  are  specified  solely  along  a  set 
of  characteristic  curves  (see  [35],  Chapter  10).  For  the  dif¬ 
ferential  equation  given  in  (4),  these  curves  are  defined  by 
all  points  {x,y)  satisfying  the  following  relation: 

Sxydy  -  ( £yy  ±  y/sfy  +  )  dx.  (6) 

To  illustrate  how  the  characteristic  curves  of  (4)  help 
formulate  elasticity  reconstruction  as  a  simplified  bound¬ 
ary  value  problem,  deformation  data  (i.e.,  all  components 
of  the  strain  tensor)  from  an  inhomogeneous  gel-based  tis¬ 
sue  equivalent  phantom  were  analyzed.  These  data  were 
collected  using  the  methods  described  in  Section  III.  The 
characteristic  curves  computed  according  to  (6)  over  a  25- 
mm  by  65-mm  area  within  the  phantom  are  presented  in 
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Fig.  1,  Characteristic  curves  of  (6)  computed  from  ultrasound  mea¬ 
surements  on  the  inhomogeneous  gel-based  phantom. 

Fig.  1.  Displayed  here  are  all  characteristics  starting  at 
the  bottom  and  left  side  of  the  region  at  equal  (approxi¬ 
mately  1.8  mm)  intervals.  If,  for  example,  the  distribution 
y)  is  given  along  parts  AB  and  AC  of  two  intersect¬ 
ing  characteristics  AB '  and  AC1 ,  then  the  reconstruction 
for  region  ABDC  reduces  to  a  classic  Goursat  boundary 
problem  for  (4).  In  contrast,  to  get  the  unique  solution  of 
(4)  within  the  region  CEF>  the  values  of  /4,  d\xjdx,  and 
dfi/dy  must  be  prescribed  along  ihe  single  line  CE  (see 
[35],  Chapter  10). 

In  this  paper  we  only  consider  reconstruction  based  on 
two  intersecting  characteristic  curves  defining  a  region  of 
interest  (ROI) ,  If  the  exact  value  of  the  elastic  modulus  is 
not  known  along  these  two  curves,  reconstruction  within 
the  ROI  is  relative.  That  is,  the  reconstructed  modulus  is 
normalized  to  the  value  along  the  intersecting  characteris¬ 
tics.  Because  the  primary  goal  of  reconstruction  is  artifact 
reduction  rather  than  exact  quantitation,  a  relative  modu¬ 
lus  image  is  sufficient.  As  discussed  in  [9],  regions  of  nearly 
constant  elastic  modulus  can  be  identified  with  edge  de¬ 
tection  operators  acting  on  strain  images. 

Using  the  characteristic  curves,  numerical  methods  can 
be  developed  to  solve  (4)  given  displacement  measure¬ 
ments.  In  practice,  however,  the  problem  of  elasticity  re¬ 
construction  is  greatly  complicated  even  in  the  linear  case 
due  to  noisy  displacement  measurements  (i.e.,  due  to  noisy 
coefficients  in  (4))  and  propagation  of  this  noise  through 
numerical  integration  within  the  ROI.  Therefore,  a  specific 
procedure  to  integrate  (4)  across  the  ROI  must  be  used  [5], 
[8]-[10],  [12],  [18]— [21]- 

To  produce  a  more  stable  reconstruction  procedure  ap¬ 
propriate  for  noisy  deformation  data,  substitute  (3)  into 
(2)  and  integrate  rather  than  differentiate.  After  eliminat¬ 


ing  the  unknown  pressure,  the  system  of  equations  reduces 
to  an  integral  equation  of  the  form: 

6{x>y)  =4 [Cyyfi  —  (^yyM)Uo  —  ieyy^)\ zq.i/d]  “b 
v 

J  [(!»),*  ~ 

VQ 

X 

-/[(7/*),y-{(7^.v}U<*x-  0,  (7) 

xo 

where  7  —  2exy  and  the  notations  f\x 0  ==  f(x o,y)  and 
f\yO  =  /(x,yo)  are  used  in  this  equation  and  below. 

The  integral  equation  in  (7)  is  expressed  as  a  functional 
5(x)y)^  as  the  goal  of  reconstruction  with  noisy  data  will 
be  to  force  6{x ,  y)  to  approach  zero  in  some  average  sense 
across  the  ROI.  In  contrast  with  (4),  this  equation  does 
not  contain  second  order  derivatives  of  the  strain.  More¬ 
over,  the  shear  strain,  and  spatial  derivatives  of  the  shear 
strain,  only  appear  in  the  integral  terms.  Because  noisy 
lateral  displacement  estimates  only  contribute  to  the  shear 
strain,  the  effects  of  measurement  error  will  be  reduced  by 
the  smoothing  action  of  the  integral  without  sacrificing 
spatial  resolution  [36].  This  type  of  processing  is  similar  to 
incompressibility  methods  in  which  noisy  lateral  displace¬ 
ment  measurements  are  smoothed  by  integrating  higher 
SNR  axial  strains  without  losing  spatial  resolution  [37], 
[38].  Therefore,  elasticity  reconstruction  by  (7)  should  be 
more  stable. 

For  large  deformations,  a  similar  functional  must  be  de¬ 
fined  from  the  general  equilibrium  equations  of  (1)  and  the 
general  Lagrangian  strain-displacement  relation.  Denoting 
the  displacement  components  u\  and  U2  by  u  and  v,  the 
unwrapped  form  of  (l)  for  the  plane  strain  state  is: 

(&xx7x  4  &xy,y)(l  4  4*  (axxU}XX  4  Pyy'tt%yy) 

4 [<Jxy,x  4  0Vy,y)w,!/  4  2a xVU,xy  5=8  0 

(0*7/, 7:  4  0ryt/,y)(I  4  W|ir)  4  {<XxxV,xx  4  &yyV,yy) 

4(ora*a;,a:  4  &xy, 4  2aXyVyXy  7=1  6* 

By  incorporating  the  stress-strain  relation  (3)  into  (8)  and 
rearranging  terms,  the  following  linear  system  of  equations 
for  derivatives  of  the  unknown  static  internal  pressure, 
and  results: 

AV(p)  =  -pB  -  F,  (9) 

where 

A(x.  y)  -  (X  +  j  «•;  J  ,  B(s,  y )  =  (bi)  =  V2t/, 

<-1,2? 

V  =  (d/dx,  d/dy),  V2  =  Q2/dx2  +  d2/dy2\ 
h  =  2[(1  4  UtX)lpi  4 

4  4  £yy^%yy  4  2 EXylL^Xy^fl\^ 

h  =  2fcXVh  +  (I  +  v%y)i>2 

4  (fisx't/xx  4  SyyV^yy  4  2 £zyV >xy)lAi 

Vh  =  {^Xx)tX  4  (flSXy)^y  1p2  —  4 
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These  formulas  contain  the  components  of  the  nonlinear 
Lagrangian  strain  tensor  [27)— [31]; 


assumption  of  incompressibility.  The  matrix  A  1  is  de¬ 
fined  as: 


E 


ij  — 


1 

2 


(10) 


(1  +  v,v  ~u,v  \ 
V  -V,x  1 
det(i4) 


(15) 


For  a  plane  strain  state,  the  strain  tensor  components  take  The  determinant  of  A  is  simply  related  to  the  metric  ten- 
the  specific  form:  sor; 


Sxx  =  U,x  +  [(u,j)J  +  (v,x)2]/2, 

evv  =  v,V  +  [(u,!/)2  +  (v,y)2  j/2,  (11) 

Exy  [w,y  H™  Vkx  H-  "t" 

If  the  magnitudes  of  the  spatial  derivatives  of  all  displace¬ 
ment  components  are  small,  the  last  (nonlinear)  term  ill 
(10)  can  be  omitted  to  produce  the  linear  strain  tensor 
of  (5). 

Again,  the  unknown  pressure  must  be  removed  from 
(9).  To  do  this,  we  solve  (9)  with  respect  to  unknowns  pfX 
and  p>y : 


dct(i4)  =  1  +  (utW  +  viV)  +  (uyXvyV  -  viXutV)  -  ^fg, 

m 

where  g  —  det(^)  i$  the  determinant  of  the  2nd  ranked 
metric  tensor  gtj. 

The  density  p  of  a  deformed  medium  is  related  to  the 
density  pQ  of  the  undeformed  medium  by  [29],  [37]: 

P  =  Po/Vo-  (17) 

For  incompressible  materials  g  =  1,  and  det(A)  ~  1. 
Therefore,  (15)  reduces  to: 


V(p)  =  mp  +  0,  (12) 

where  a{x,y)  =  (a,:)  =  — 0{x,y)  =  (A)  =  -A-'F, 
i  —  1,2.  Then,  integrating  the  first  of  these  equations  along 
x  and  the  second  along  two  simultaneous  expressions  for 
the  unknown  pressure  are  produced: 


A~1  _  A  +  V.V  ~u,v  ^ 

\  -V.x  1  +  u,xj 


(18) 


Substituting  (18)  into  (14)  and  integrating,  we  obtain 
an  expression  for  the  reconstruction  functional  in  the  non¬ 
linear  case,  as  seen  in  (19)  (top  of  next  page),  where 


p(x,y)  =  w 


p{x,y)  =  ip‘2 


T. 

p{x  o.  u)  +  J 


(h 

<pi 


dx 


V 


P(*>  Vo)  +  J 


wo 


02 

<Pi 


dy 


(13.1) 


(13.2) 


where  <pi(x,y)  =  exp  j/“o  q,  (x,  j/)dx|  and  <Pi(.x,y)  = 

“p{/:  a-i(z,  y)dy}  do  not.  depend  on  n(x,  y).  The  term 
p(XiVo)  in  (13.2)  can  be  obtained  from  (13,1).  Similarly, 
the  p{xo%y)  in  (13.1)  can  be  obtained  from  (13.2).  Conse¬ 
quently,  the  set  of  simultaneous  equations  describing  the 
unknown  pressure  reduces  to  a  single  equation, 


'fix  f  —dx  -ifi-2  (pi  I  —dx  -  <p2  [  —dy 

J  V 1  J  <Pl  J  <J>2 

x0  L  Vo  J  UQ 


+  V  i 


(f2 

L  Vo 


vo  vq 

-  Gpo  =  0,  (14) 


where  G(x,y)  ~  [(p2(^y)fpi(x7y0)  -  tp\{x,y)y2(xo,y]\ 
does  not  depend  on  p(x,  y),  and  p0  =  p(x 0,  yQ). 

The  expression  given  in  (14)  can  be  used  to  reconstruct 
the  shear  modulus  given  large  deformations,  Before  defin¬ 
ing  a  functional  similar  to  (7)  for  the  nonlinear  case,  we 
consider  a  simplification  of  the  matrix  A”1  based  on  an 


€  —  Eyy  txx,  7  —  2S 

9\  =  2{ekvv(l  +  vM)  -  viVVuiV) 

+  7k«v(l  +  Vty)  -  V,=yUy}, 

92  =  “b  'U,x)  t^JTSC^.a?] 

*b  "h  ^,a:)  “ 

In  the  limit  of  small  displacements,  a  =?  0,  =  1,  i  =  1, 2, 

G(^7^/)  =  0,  and  all  second  order  terms  in  p*  can  be  ig¬ 
nored.  For  this  case,  (19)  reduces  to  (7).  Consequently,  (19) 
can  be  used  as  the  reconstruction  functional  over  a  wide 
range  of  internal  deformations.  In  contrast  to  (4),  the  dif¬ 
ferential  equation  for  the  nonlinear  case  after  eliminating  p 
contains  4th  order  displacement  derivatives,  Therefore,  an 
integral  representation  for  the  nonlinear  case  is  also  much 
more  appropriate  given  realistic  measurement  noise. 

To  evaluate  either  (7)  or  (19)  7  spatial  derivatives  of 
both  the  axial  and  lateral  displacement  components  are 
needed.  It  is  well  known,  however,  that  lateral  displace¬ 
ment  measurements  are  much  noisier  than  axial  ones  if 
ultrasound  is  used  to  track  internal  deformations  [7],  [37], 
[38] .  In  addition  to  smoothing  from  the  integral  represen¬ 
tation  of  the  equilibrium  equation,  incompressibility  pro¬ 
cessing  Can  be  used  in  both  linear  and  nonlinear  cases  to 
further  reduce  the  influence  of  noisy  displacement  mea¬ 
surements  [37],  [38],  All  results  presented  in  Section  III 
used  lateral  displacement  estimates  obtained  from  incom¬ 
pressibility  processing. 

To  solve  for  the  unknown  shear  elasticity  distribu¬ 
tion  y)  given  all  measured  displacements  and  strains 
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y)  —  2  je/i  [ep  +  {£xxV)\xoVl]\yQV2  -  [ep  -  {^yyP)\yQ^P2]\XQ^Pl^  + 

[(7m),x  +  \ig^P2l(^y  |  y>2  -  <  V2  J [(7a0,®  +  M2\<P2ldy 

)  V  2/0 

J[(nv),v  +  Wilvr1*}  fi  -  |<pi  J [(7 m),v  + 


<P1  - 


<^2  -  Gpo  =  o,  (19) 


y  o 


within  the  ROI,  a  global  minimization  procedure  was  used 
[5],  [8]— [10] ,  [12],  [18]— [21].  This  numerical  technique  is 
general  and  applies  both  to  linear  (7)  and  nonlinear  (19) 
functionals.  For  a  given  distribution  p(x,  y)  along  the  ROI 
boundaries,  the  error  functional  5(x,  y)  must  be  minimized 
across  the  ROI  in  some  general  way.  Here,  the  specific  dis¬ 
tribution  {iij  =  fi(xi,yj)  is  sought  which  minimizes  the 
total  error: 


where  (£;,%•),  0  <  i  <  iV,  0  <  j  <  M  is  a  rectangular 
grid  covering  the  ROI.  The  integral  is  approximated  by 
summing  6(xi,yi)  over  all  grid  points. 

In  discrete  form,  when  differential  and  integral  opera¬ 
tors  are  replaced  by  discrete-space  equivalents,  simultane¬ 
ous  minimization  of  D  with  respect  to  all  undetermined 
fJLij  yields  a  set  of  linear  algebraic  equations.  In  theory, 
therefore,  reconstruction  is  primarily  a  matrix  inversion. 
In  practice,  however,  the  matrix  is  poorly  conditioned;  it 
is  very  difficult  to  produce  a  stable  inverse  using  noisy,  ex¬ 
perimental  data.  Alternatively,  error  minimization  can  be 
performed  with  an  iterative  approach,  as  discussed  below. 

For  the  nonlinear  case,  the  unknown  scalar  value  po  — 
p{x0,yo)  in  (19)  must  be  estimated  prior  to  general  recon¬ 
struction.  This  term  is  estimated  independently  by  mini¬ 
mizing  the  error: 


within  a  thin  region  including  the  ROI  boundaries,  x  = 
x0  and  y  =  y0 ,  S*  =  (x0  <  x  <  xi ,y0  <  y  <  yM)  U 
{%o  <  x  <  XN,yo  <  y  <  yi),  in  which  the  distribution  of 
p(x,  y)  is  assumed  known  given  the  specified  values  along 
the  boundaries  themselves. 

The  specific  iterative  procedure  used  here  to  compute 
the  distribution  pij  =  p(xi,yj)  minimizing  either  (7)  or 
(19)  over  the  ROI  is  based  on  a  gradient  method 

/4+1  =  4  -  XijdD/dmj,  (22) 

where  k  is  the  iteration  index.  It  starts  with  a  trial  so¬ 
lution  Then,  the  error  D  is  minimized  by  varying 


the  unknown  shear  modulus  /x  at  only  one  given  grid  point 
(ij).  This  procedure  is  repeated  for  each  (i,  j)  grid  point 
([35],  Chapter  20).  The  iterative  parameters  A £■,  which  de¬ 
termine  the  step  size  of  the  gradient  method,  were  chosen 
based  on  three  estimates  of  D.  That  is,  the  minimum  of 
D  was  locally  predicted  using  a  second  order  polynomial 
approximation  of  D  as  a  function  of  at  each  grid  point 
under  the  restriction  of  a  decreasing  error  [5],  [8]— [10],  [12], 
[18]— [21] ,  [35].  Then,  a  global  linear  predictor  was  used 
to  update  all  simultaneously.  This  reduced  the  oscil¬ 
latory  nature  of  convergence.  If  the  total  error  remained 
nearly  constant  at  a  given  step  k  based  on  the  global  lin¬ 
ear  predictor  values,  then  the  iterative  step  sizes  A ^  were 
again  selected  separately  using  the  local  quadratic  predic¬ 
tor  as  described  above.  By  “ping-ponging”  between  local 
and  global  criteria  in  this  way,  large  oscillations  as  a  func¬ 
tion  of  iteration  index  were  greatly  reduced,  thus  speeding 
convergence. 

All  spatial  derivatives  in  the  reconstruction  equations 
were  replaced  with  2nd  order  finite  differences  over  the 
same  grid  (a^,  yj).  Because  2nd  order  finite  differences  use 
only  information  from  neighboring  pixels,  error  computa¬ 
tions  were  optimized  so  that  for  each  pioj0  update,  the 
error  Sij  =  5(xi,yj)  was  computed  only  for  \i  -  z0|  <  1  or 
|  j  -  jo  I  <  1-  Computations  continued  until  the  total  error 
reached  a  stable  plateau  \Dk +1  —  Dk\/Dk  <  So-  For  all 
results  presented  below,  a  value  of  So  =  10“6  was  used, 
and  a  homogeneous  medium  was  the  trial  solution.  It  was 
observed  that,  with  different  displacement  fields  and  differ¬ 
ent  grid  sizes,  the  algorithm  demonstrates  approximately 
exponential  convergence.  By  decreasing  the  grid  size,  the 
rate  of  convergence  decreased  almost  linearly.  Typical  re¬ 
construction  times  on  a  low-end  SPARC  2  workstation  for 
the  results  presented  below  were  less  than  5  minutes  for 
the  nonlinear  functional  (19),  and  even  less  for  the  linear 
functional  (7).  These  algorithms,  however,  can  be  further 
optimized  for  time  performance,  but  this  is  beyond  the 
scope  of  this  paper. 


III.  Experimental  Methods  and  Materials 

All  measurements  were  made  with  an  ultrasound-based 
deformational  imaging  system  similar  to  the  one  presented 
in  [24].  A  38-mm  wide,  linear  array  transducer  operating  at 
5  MHz  was  used  for  all  studies.  The  array  was  driven  with 
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an  Ultramark-9  (ATL  Corp,  Bothell,  WA)  real-time  ultra¬ 
sound  system  operating  in  conventional  B-scan  mode.  The 
digital  radio  frequency  (RF)  signal  output  by  the  beam- 
former  was  captured  before  subsequent  processing  and  dis¬ 
play  by  the  Ultramark-9  back-end.  By  buffering  RF  data  at 
the  beamformer  output  for  data  capture  with  an  external 
device,  live  images  could  be  viewed  during  data  capture. 
About  120  consecutive  frames  of  real-time  RF  data  were 
stored  using  a  digital  data  capture  system  constructed  in 
our  lab.  At  a  typical  35  Hz  frame  rate,  this  represents  al¬ 
most  4  seconds  of  phase  sensitive  ultrasound  data  that  can 
be  used  for  sensitive  speckle  tracking. 

Measurements  were  made  on  two  different  gel-based 
phantoms.  These  phantoms  were  constructed  using  the 
procedures  described  in  [5]— [10].  Both  phantoms  were 
100  mm  wide  by  140  mm  long,  where  all  applied  surface 
deformations  were  vertical,  and  the  depth  direction  in  the 
ultrasound  imaging  plane  also  was  vertical.  The  first  phan¬ 
tom  was  homogeneous  and  measured  117  mm  in  height. 
The  second  was  80  mm  high  and  included  a  single  cylin¬ 
drical  hard  inclusion  near  the  bottom  of  the  phantom.  This 
inclusion  in  the  inhomogeneous  phantom  was  18  mm  in  di¬ 
ameter  and  was  oriented  so  that  the  longitudinal  axis  of  the 
cylinder  was  perpendicular  to  both  the  ultrasound  image 
plane  and  the  direction  of  the  applied  surface  deformation. 
The  inclusion  was  constructed  from  a  higher  concentration 
§cl*  This  gel  had  a  shear  modulus  about  2.5  times  larger 
than  the  surrounding  material  as  estimated  from  indepen¬ 
dent  measurements  of  the  elastic  modulus  using  the  system 
described  in  [25]. 


Fig.  2.  Measured  axial  displacement  images  from  the  homogeneous 
phantom  (left)  and  phantom  with  single  hard  inclusion  at  the  bot¬ 
tom  (right). 


Both  phantoms  were  vertically  deformed  with  a 
12.5  mm  thick  Plexiglas  plate  attached  to  a  manually  con¬ 
trolled,  one-dimensional  motion  axis.  The  plate  measured 
125  mm  by  70  mm  in  cross  section,  and  almost  completely 
covered  the  top  surface  of  either  phantom.  Such  a  large 
plate  ensured  that  a  plane  strain  state  was  reasonably  ap¬ 
proximated  for  both  phantoms.  A  hole  was  cut  into  the 
center  to  mount  the  imaging  array.  Once  the  array  was 
properly  secured,  the  bottom  surface  of  the  plate  main¬ 
tained  continuous  contact  with  the  top  surface  of  the  phan¬ 
tom.  Deformations  were  applied  by  smoothly  turning  the 
gear  to  move  the  plate  vertically  over  a  distance  up  to 
30  mm  during  the  4  second  data  capture  period.  Conse¬ 
quently,  very  large  vertical  displacements  could  be  applied 
during  continuous  ultrasound  data  capture. 

To  test  both  linear  and  nonlinear  reconstruction  pro¬ 
cedures  over  a  wide  deformational  range,  two  data  sets 
were  recorded  for  each  phantom.  Internal  displacements 
and  strains  were  imaged  in  the  homogeneous  phantom  for 
applied  surface  displacements  of  0.7  mm  (0.6%  average  ver¬ 
tical  strain),  representing  a  small  deformation,  and  7.6  mm 
(6.5%  average  vertical  strain),  representing  a  fairly  large 
deformation.  For  the  phantom  with  a  single  hard  inclu¬ 
sion,  larger  surface  displacements  were  used  to  produce 
strains  within  the  hard  inclusion  comparable  to  the  av¬ 
erage  strains  in  the  homogeneous  phantom.  The  two  data 
sets  recorded  on  this  phantom  used  a  surface  displacement 
of  2.7  mm  (3.4%  average  vertical  strain)  for  the  small  de¬ 


formation  case,  and  12.8  mm  (16%  average  vertical  strain) 
for  the  large  deformation  case. 

All  displacement  and  strain  images  were  computed  from 
RF  ultrasound  data  using  the  speckle  tracking  procedures 
described  in  [39],  [40].  Based  on  spatial  autocorrelation 
analysis  of  the  axial,  vertical  strain  image  in  the  homo¬ 
geneous  phantom  [39],  the  spatial  resolution  of  these  im¬ 
ages  was  estimated  to  be  about  1.8  mm.  Consequently, 
the  grid  used  for  all  reconstructions  had  equal  1.8  mm 
spacing  in  both  directions  (i.e.,  Ax  =  Ay  =  1.8  mm). 
The  shear  modulus  was  reconstructed  in  all  cases  within 
a  25.1  mm  x  66.4  mm  rectangular  ROI  located  near  the 
vertical  center  line  of  the  ultrasound  image. 

For  ultrasound  speckle  tracking  using  RF  data,  lateral 
displacement  estimates  exhibit  significantly  lower  SNR 
than  axial  (vertical)  estimates  [38].  To  overcome  this  limi¬ 
tation,  incompressibility  processing  methods  have  been  de¬ 
veloped  for  linear  and  nonlinear  cases  [37],  [38].  For  large 
deformations,  the  incompressibility  condition  is: 

(1  +  v,y)u, x  -  v,xu,y  +  v.y  =  0.  (23) 

As  demonstrated  in  [37],  the  measured  lateral  displace¬ 
ment  image  uexp(x,  y)  is  used  as  a  reference  during  lateral 
displacement  computation.  For  the  images  presented  here, 
the  reference  area  was  a  rectangle  10.5-mm  wide  located 
near  the  vertical  central  line  of  the  ultrasound  image.  In 
contrast  with  [37],  where  a  polynomial  approximation  of 
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Fig.  3-  Measured  lateral  displacements  images  for  the  homogeneous  phantom  (far  left)  and  the  phantom  with  a  single  hard  inclusion  (second 
from  left),  Incompressibility  processed  lateral  displacement  images  for  the  homogeneous  phantom  (second  from  right)  and  the  phantom  with 
a  single  hard  inclusion  (far  right) . 


the  unknown  lateral  displacement  along  an  entire  single 
line  yielded  a  solution  of  (23),  a  more  flexible  approach 
was  used  here;  (23)  was  solved  for  a  set  of  particular  solu¬ 
tions  of  the  form  up(x 0,  y )  =  yv,p~  1, 2, 3  (i.e,,  3rd  order) 
within  the  region  j/o  -  A  <  y  <  j/0  -h  A  for  every  interior 
point  (xotVq)-  Therefore,  a  set  of  particular  solutions  of 
(23)  was  obtained  within  this  region,  in  which  the  lateral 
displacement  u(x,  y)  is  a  linear  combination  of  these  par¬ 
ticular  solutions.  The  three  unknown  coefficients  defining 
the  linear  combination  of  particular  solutions  were  found 
by  minimizing  the  total  error : 

(tt  -  u*xp)2ds  (24) 

Sn 

across  the  corresponding  area  =  (x3  <  x  <  Xf , 

2/o  —  A  <  y  <  yo  +  A)  of  the  referenced  region. 

Prior  to  reconstruction,  both  measured  axial  displace¬ 
ments  and  incompressibility  processed  lateral  displace¬ 
ments  were  filtered  with  a  two-dimensional  Hamming  func¬ 
tion,  further  reducing  the  spatial  resolution  of  displace¬ 
ment  images  and  subsequent  elasticity  reconstructions  to 
about  2.5  mm. 

IV.  Results 

Measured  vertical  (axial)  displacement  images  over  the 
ROI  for  the  large  deformation  case  are  presented  in  Fig.  2 


for  the  homogeneous  phantom  [Fig.  2(a)]  and  the  inho¬ 
mogeneous  phantom  [Fig.  2(b)].  In  all  these  images  the 
transducer,  and  hence  the  reference  for  all  displacement 
measurements,  is  at  the  top.  Consequently  vertical  dis¬ 
placements  arc  zero  at  the  top,  and  motion  is  toward  the 
transducer  (i.e.,  negative  vertical  motion).  The  display  dy¬ 
namic  range  for  Fig.  2(a)  is  —5.6  mm  to  —0.4  mm,  and  for 
Fig.  2(b)  it  is  -12,6  mm  to  —1-0  mm.  Images  of  the  mea¬ 
sured  lateral  displacement  are  presented  in  Fig.  3  under  the 
same  conditions  for  the  homogeneous  phantom  [Fig,  3(a)] 
and  the  inhomogeneous  phantom  [Fig.  3(b)].  The  display 
dynamic  range  for  the  homogeneous  phantom  is  —  1  mm 
to  1  mm,  and  for  the  inhomogeneous  phantom  is  -1  mm 
to  4  mm,  in  which  black  represents  motion  to  the  left  and 
white  is  motion  to  the  right-  (The  bottom  of  the  phantom 
shifted  to  the  right  for  the  inhomogeneous  case,  and  con¬ 
sequently  an  asymmetric  display  dynamic  range  is  used.) 
Using  the  measured  vertical  displacements  of  Fig,  2,  the 
same  lateral  displacement  images  after  incompressibility 
processing  are  presented  in  Fig,  3(c)  for  the  homogeneous 
phantom,  and  Fig,  3(d)  for  the  inhomogeneous  phantom. 
Images  of  the  shear  modulus  reconstructed  using  (7) 
and  (19)  are  presented  in  Fig.  4  for  the  homogeneous 
phantom.  Fig.  4(a)  represents  the  linear  reconstruction  for 
the  small  deformation  case,  whereas  Fig,  4(b)  represents 
the  nonlinear  reconstruction.  Reconstructions  for  the  large 
deformation  case  are  presented  in  the  lower  two  panels, 
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where  Fig.  4(c)  is  the  linear  reconstruction  and  Fig.  4( 
is  the  nonlinear  one.  These  relative  modulus  images  a 
presented  over  a  logarithmic  gray  scale  (0.47  to  2  12) 
illustrated  on  the  right.  The  mean  value  of  the  norm 
ized  (i.e  normahzed  to  the  average  elastic  modulus  alo 
the  boundaries  of  the  ROI)  reconstructed  elasticity  d 
tribution  within  the  ROI  is  1.032  for  Fig.  4(a)  1.022  f 
Fig.  4(c) :  0.991  for  Fig.  4(b)  and  0.996  for  Fig  4(d)  T1 
standard  deviation  is  0.048  and  0.023  for  Figs.  4(a)  and  (c 
and  0.033,  and  0.019  for  Figs.  4(b)  and  (d),  respective! 
Artifacts  in  the  nonlinear  reconstructions  are  no  great< 
than  those  in  the  linear  reconstructions  even  for  small  d< 
formations.  Horizontal  artifacts,  created  by  low  SNR  lai 
era  displacement  measurements,  are  present  in  both  lir 
ear  and  nonlinear  reconstructions.  With  depth,  these  ai 
tifacts  have  more  energy  at  lower  spatial  frequencies  (i  e 
are  br  oader  at  the  bottom  of  the  image)  as  the  lateral  poin 
spread  function  of  the  ultrasound  system  broadens. 

Images  of  the  shear  modulus  reconstructed  using  (7 
and  19)  are  presented  in  Fig.  5  for  the  inhomogeneom 
phantom  Fig.  5(a)  represents  the  linear  reconstruction  foi 
the  small  deformation  case,  whereas  Fig.  5(b)  represents 
the  nonlinear  reconstruction.  Reconstructions  for  the  large 
deformation  case  are  presented  in  the  lower  two  panels, 
where  Fig.  5(c)  is  the  linear  reconstruction  and  Fig.  5(d) 
1  6  nc™^near  one-  Exactly  the  same  logarithmic  gray 
scale  as  Fig.  4  was  used  here.  Again,  horizontal  artifacts 


broaden  with  depth,  becoming  especially  noticeable  at  the 
oot tom  or  these  images. 

In  all  reconstructions,  the  uniform  mechanical  bound- 

the  ROT  R°n  M  =  1  WaS  358111116(1  aIon§  the  edge  of 
the  ROI.  Because  the  equilibrium  equation  is  hyperbolic 

any  errors  in  this  assumption  produce  noise  across  the  en¬ 
tire  reconstruction  where  the  relative  magnitude  of  this 
se  (i.e.,  | n  n\/n  with  jj.  the  reconstructed  image  and 

M  the  true  image)  is  about  the  same  as  the  magnitude  of 
the  boundary  condition  error.  This  analysis  was  confirmed 
y  simulating  numerous  elasticity  reconstructions  of  a  sin¬ 
gle  hard  inclusion  in  an  otherwise  homogeneous  medium, 
where  boundary  errors  were  modeled  by  a  10-term  Fourier 
series  with  random  coefficients. 


The  results  presented  in  Figs.  4  and  5  demonstrate  that 
6  quality  of  the  image  reconstructed  with  the  i!geomet- 
nc  nonlinear  model  rivals  that  of  the  reconstruction  based 
on  a  linear  model,  even  though  high  order  displacement 
derivatives  are  used.. Moreover,  nonlinear  processing  pro- 
vi  es  more  consistent  elasticity  reconstructions  even  if  hor¬ 
izontal  streak  artifacts  are  slightly  elevated  compared  to 
linear  processing.  To  illustrate  this,  the  elasticity  distri- 
u  ion  reconstructed  along  the  central  vertical  line  of  the 
inhomogeneous  phantom  is  presented  in  Fig.  6  for  linear 
Eig.  6(a)]  and  nonlinear  [Fig.  6(b)]  models.  Clearly,  non¬ 
linear  processing  produces  more  consistent  reconstructions 
independent  of  the  magnitude  of  the  external  deforma- 
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Fig.  5.  Elasticity  distribution  for  the  phantom  with  a  cylindrical  hard  inclusion  at  the  bottom  reconstructed  by  linear  processing  (far  left)  and 
(second  from  right),  and  nonlinear  processing  (second  from  left)  and  (far  right).  Displacement  images  from  a  3.49c  mean  vertical  deformation 
were  used  in  (far  left)  and  (second  from  left).  Displacement  images  from  a  16%  mean  vertical  deformation  were  used  in  (second  from  ri<dit) 
and  (far  right).  G  ' 


tion.  This  property  holds  not  only  for  single  line  profiles. 
For  the  homogeneous  phantom,  the  mean  squared  differ¬ 
ence  between  small  and  large  deformation  images  com¬ 
puted  over  the  entire  ROI  is  0.115  for  the  linear  recon¬ 
struction  [Figs.  4(a)  and  (c)]  and  0.076  for  the  nonlinear 
reconstruction  [Figs.  4(b)  and  (d)].  Similarly,  for  the  inclu¬ 
sion  phantom,  the  mean  squared  difference  between  small 
and  large  deformation  images  computed  over  the  top  half 
of  the  ROI  is  0.086,  and  over  the  bottom  half  of  the  phan¬ 
tom  it  is  0.181  for  the  linear  reconstruction  [Figs.  5(a) 
and  (c)],  and  0.071  for  the  top  half  of  the  phantom  and 
0.152  for  the  bottom  half  of  the  nonlinear  reconstruction 
[Figs.  5(b)  and  (d)]. 


V.  Discussion 

The  results  of  elasticity  reconstructions  presented  above 
show  that  a  nonlinear  (geometric)  model  can  be  used  to 
reconstruct  the  shear  (or  Young’s)  elastic  modulus  based 
on  internal  displacement  and  strain  fields  computed  from 
real-time  ultrasound  data.  Despite  the  high  order  spatial 
derivatives  of  both  displacement  components  required  for 
this  processing,  the  image  quality  of  reconstructions  based 
on  a  nonlinear  model  rival  those  based  on  a  simpler,  lin¬ 
ear  model.  Moreover,  nonlinear  processing  provides  more 


consistent  results  so  that  elasticity  reconstruction  can  be 
performed  over  a  wide  range  of  external  loading. 

The  specific  numerical  methods  developed  here  exploit 
an  integral  rather  than  differential  representation  of  the 
elasticity  equations.  Given  noisy  displacement  and  strain 
estimates,  this  approach  is  more  robust  and  leads  to  more 
stable  reconstructions.  Moreover,  coupled  with  an  iterative 
procedure  using  both  local  and  global  error  predictors  to 
accelerate  convergence,  it  can  produce  shear  modulus  im¬ 
ages  within  a  few  minutes  on  a  low-end,  general  purpose 
computer.  Future  work  will  optimize  software  running  on 
more  powerful  workstations  to  reduce  reconstruction  times 
to  a  few  tens  of  seconds.  Such  times  are  appropriate  for 
clinical  applications  using  real-time  ultrasound  data  cap¬ 
ture  and  specially  constructed  hardware  for  speckle  track¬ 
ing. 

The  fundamental  hypothesis  of  this  study  is  that  re¬ 
construction  can  greatly  reduce  artifacts  in  displacement 
and  strain  images  due  to  global  boundary  conditions.  To 
illustrate  this  point  for  the  specific  example  discussed  here, 
the  results  of  a  simple  1-D  elasticity  reconstruction  are  pre¬ 
sented  in  Fig.  7.  Fig.  7  shows  the  normalized  distribution  of 
l/viV  for  the  homogeneous  phantom  [Fig.  7(a)],  and  the  in¬ 
homogeneous  phantom  [Fig.  7(b)]  over  precisely  the  same 
display  dynamic  range  used  in  Figs.  4  and  5.  Because  the 
homogeneous  phantom  was  constructed  to  yield  a  nearly 
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Fig.  7.  Normalized  distribution  of  l/vty  for  the  homogeneous  phan¬ 
tom  (left),  and  the  phantom  with  a  single  hard  inclusion  (right). 


Fig.  6.  Elasticity  distribution  along  central  vertical  line  of  the  inho¬ 
mogeneous  phantom  images,  reconstructed  by  linear  (left)  and  non¬ 
linear  (right)  processing. 


uniform  strain  distribution  over  the  ROI,  the  1-D  recon¬ 
struction  is  almost  perfect.  That  is,  measurement  noise  is 
not  amplified  for  this  type  of  processing,  leading  to  a  clean 
reconstruction.  In  contrast,  this  simple  approach  produces 
significant  artifacts  in  images  of  an  inhomogeneous  object. 
Clearly,  as  the  object  becomes  more  complex,  more  accu¬ 
rate  reconstruction  is  required  even  when  uniform,  one¬ 
dimensional  loading  is  used. 

All  the  results  presented  in  this  paper  were  computed 
assuming  a  linear  stress-strain  relation  (3).  This  assump¬ 
tion  is  almost  perfect  over  a  wide  deformation  range  for  the 
gel-based,  tissue-like  phantoms  used  here  [41].  Real  tissue, 
however,  exhibits  nonlinear  behavior  (i.e.,  material  nonlin¬ 
earity)  even  for  simple  types  of  external  loading  if  the  de¬ 
formation  is  significant  [25],  [34],  [42].  Because  ultrasound 
data  can  be  acquired  over  a  wide  deformation  range,  there 
is  the  possibility  that  the  shear  modulus  can  be  recon¬ 
structed  at  different  strain  magnitudes.  Over  a  limited  de¬ 
formation  range,  the  stress-strain  curve  can  be  considered 
linear,  but  with  an  elastic  modulus  that  depends  on  the 
instantaneous  strain  magnitude.  Displacement  and  strain 
data  for  an  ex  vivo  model  of  kidney  transplant  rejection 
acquired  over  a  large  deformation  range  have  been  ana¬ 
lyzed  with  a  piecewise  linear  approximation  to  produce 


images  related  to  the  material  nonlinearity  of  the  kidney 
[24].  Future  work  will  combine  such  measurements  with 
the  methods  developed  here  to  image  the  nonlinear  elastic 
properties  of  tissue. 
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Appendix 


Assume  that  two  displacement  components  are  mea¬ 
sured  with  high  precision  within  the  volume;  also,  as¬ 
sume  that  the  third  component  is  estimated  less  accu¬ 
rately,  as,  for  example,  in  elasticity  measurements  using 

MRI  [18],  [21]. 

For  a  general  three-dimensional  strain  state  described 
by  (1)  and  (3),  the  following  system  of  equations  analogous 
to  (9)  is  produced: 


AV(p)  =  -pB  -  F, 


(Al! 
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where 


(l  +  wl,l  ul,2  ^1,3  \ 

u2,l  1  +  1*2,2  1*2,3  I  , 

w3,l  **3,2  1+1*3, 3/ 

B{x,  y)  =  (bi)  =  V2U, 

F  =  ( fi ) 


^  =  (ip i)  =  (p£li),l  +  (/U£2i),2  +  (/^3*),3,  *  =  1, 2, 3; 

V  =  (d/dx1,d/dx2,d/dx3), 
v2  =  d2/dx2  +  d2/dx2  +  d2/dx2. 

Note  that  the  nonlinear  form  of  the  strain  tensor  (10)  must 
be  used. 

To  solve  (Al)  with  respect  to  unknowns  dp/d: ri, 
dp/dx 2,  and  dp/dx 3,  we  first  compute  the  determinant 
of  matrix  A,  which  is 


where  g  —  det(5y)  is  the  determinant  of  the  2nd  ranked 
metric  tensor  gij. 

Note  that  for  incompressible  materials  g  =  1,  and, 
det(A)  =  1  which  greatly  simplifies  the  inversion  of  (Al): 

V  (p)  =  ap  +  0,  (A2) 

where  a{x,y)  =  (<*)  =  -A~1B,  f3{x,y )  =  (ft)  =  -A~lF, 
i  =  1,2,3. 

By  integrating  each  component  equation  of  (A2)  along 
its  corresponding  coordinate,  we  obtain: 

p(xi,x2,x3)  =  <fiip(xi,x2,x3)  +  Fi 

=  V2P{xi ,  x%,  x3)  +  F2  (A3) 
=  V3P{xi,x2,x3)  +  F3 


\  Aty  p 

3  3 

+  £  (1  - 

{ 

i=\  i,j~l 

where 


<Pi(xi,X21X3)  =  exp 


Xi 

J  ai(xi1X2,x3)dxi 


Fi(x  i,x2,x3)  =  ^J~dxiJ  <pu  i  =  1,2,3. 

Using  a  procedure  paralleling  that  used  to  produce  (14), 
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the  three-dimensional  case  reduces  ter 

* 

+  VW2|xo)p0  +  (if2F3  |xo  +  <p3F2\xo) 

+  (F2  +  F3)]|x?^  +  2  Fi  = 

+  FsFi  Ixg  )Pq  +  (<FiF3|xo  +  v?3Fi|xo) 

+  (Fi  +  F3)\\x%V2  +  2  F2  = 

[(V’l^lx?  +  P2<Fl|x°)PO  +  (<FlF2|xo  +v?2Ti|xo) 

+  (Fi  +  F2)]|xo^3  +  2F3, 

where  p0  =  p(x°,x%,x3).  For  a  plane  strain  state  (A4) 
reduces  to  (14). 

Again,  in  the  limit  of  small  displacements,  a  =  0,  <pi  = 
1,  z  =  1, 2, 3,  (A4)  can  be  greatly  simplified: 

(FsIxO  +  F2|X0  +F2  +  F3)|x0  +  2F,  = 

(F3|xo  +  Fi|xo  +FX+  F3)|xc  +  2  F2  = 

(F2|xo  +  Fi|xo  +  Fi  +  F2)|xo  +  2F3,  (A5) 

where 


Fi  =  2  < 


F>  =  2  l 


F3  =  2  < 


(m£ii)Ix“  _  P-£ li  -  J  [(pe  12), 2  -  (ja£i3),3]dxi 

x° 
x2 

(^22)|xo  -  F^22  -  J[{pe  12), 1  -  0u£23),3]d:T2  1  , 

2*0 

*3 

(M£33)|x»  -  P£33  ~  J [(/t£l3),l  -  (^23), 2 }dx3 


Equations  (A4)  and  (A5)  do  not  contain  high  order  spatial 
displacement  derivatives,  compared  with  their  equivalent 
differential  equations,  and  therefore,  elasticity  reconstruc¬ 
tion  by  (A4)  and  (A5)  should  be  more  stable.  Again,  (A4) 
and  (A5)  show  that  spatial  derivatives  of  all  displacement 
components  are  needed  in  general  for  elasticity  reconstruc¬ 
tion  in  both  linear  and  nonlinear  cases,  but  any  one  of 
the  three  displacement  components  can  be  reconstructed 
using  incompressibility  processing  based  on  the  relation 
(det(A)  =  1)  [37]. 
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